#include <BCIntegrate.h>
Inherits BCEngineMCMC.
Inherited by BCModel.
Inheritance diagram for BCIntegrate:
Member functions (miscellaneous methods) | |
| void | DeleteVarList () |
| void | ResetVarlist (int v) |
| void | UnsetVar (int index) |
| virtual double | Eval (std::vector< double > x) |
| virtual double | LogEval (std::vector< double > x) |
| virtual double | EvalSampling (std::vector< double > x) |
| double | LogEvalSampling (std::vector< double > x) |
| double | EvalPrint (std::vector< double > x) |
| virtual double | FitFunction (std::vector< double > x, std::vector< double > parameters) |
| double | Integrate () |
| double | IntegralMC (std::vector< double > x, int *varlist) |
| double | IntegralMC (std::vector< double > x) |
| double | IntegralMetro (std::vector< double > x) |
| double | IntegralImportance (std::vector< double > x) |
| double | CubaIntegrate (int method, std::vector< double > parameters_double, std::vector< double > parameters_int) |
| double | CubaIntegrate () |
| TH1D * | Marginalize (BCParameter *parameter) |
| TH2D * | Marginalize (BCParameter *parameter1, BCParameter *parameter2) |
| TH1D * | MarginalizeByIntegrate (BCParameter *parameter) |
| TH2D * | MarginalizeByIntegrate (BCParameter *parameter1, BCParameter *parameter2) |
| TH1D * | MarginalizeByMetro (BCParameter *parameter) |
| TH2D * | MarginalizeByMetro (BCParameter *parameter1, BCParameter *parameter2) |
| int | MarginalizeAllByMetro (const char *name) |
| TH1D * | GetH1D (int parIndex) |
| int | GetH2DIndex (int parIndex1, int parIndex2) |
| TH2D * | GetH2D (int parIndex1, int parIndex2) |
| void | InitMetro () |
| void | SAInitialize () |
| virtual void | FindModeMinuit (std::vector< double > start=std::vector< double >(0), int printlevel=1) |
| void | FindModeMCMC () |
| void | FindModeSA (std::vector< double > start=std::vector< double >(0)) |
| double | SATemperature (double t) |
| double | SATemperatureBoltzmann (double t) |
| double | SATemperatureCauchy (double t) |
| virtual double | SATemperatureCustom (double t) |
| std::vector< double > | GetProposalPointSA (std::vector< double > x, int t) |
| std::vector< double > | GetProposalPointSABoltzmann (std::vector< double > x, int t) |
| std::vector< double > | GetProposalPointSACauchy (std::vector< double > x, int t) |
| virtual std::vector< double > | GetProposalPointSACustom (std::vector< double > x, int t) |
| std::vector< double > | SAHelperGetRandomPointOnHypersphere () |
| double | SAHelperGetRadialCauchy () |
| double | SAHelperSinusToNIntegral (int dim, double theta) |
| virtual void | MCMCUserIterationInterface () |
| static void | CubaIntegrand (const int *ndim, const double xx[], const int *ncomp, double ff[]) |
| static void | FCNLikelihood (int &npar, double *grad, double &fval, double *par, int flag) |
Public Types | |
Enumerators | |
| enum | BCIntegrationMethod { kIntMonteCarlo, kIntImportance, kIntMetropolis, kIntCuba } |
| enum | BCMarginalizationMethod { kMargMonteCarlo, kMargMetropolis } |
| enum | BCOptimizationMethod { kOptSA, kOptMetropolis, kOptMinuit } |
| enum | BCSASchedule { kSACauchy, kSABoltzmann, kSACustom } |
Public Member Functions | |
Member functions (get) | |
| BCIntegrate::BCIntegrationMethod | GetIntegrationMethod () |
| BCIntegrate::BCMarginalizationMethod | GetMarginalizationMethod () |
| BCIntegrate::BCOptimizationMethod | GetOptimizationMethod () |
| BCIntegrate::BCOptimizationMethod | GetOptimizationMethodMode () |
| BCIntegrate::BCSASchedule | GetSASchedule () |
| void | GetRandomVector (std::vector< double > &x) |
| virtual void | GetRandomVectorMetro (std::vector< double > &x) |
| double | GetRandomPoint (std::vector< double > &x) |
| double | GetRandomPointImportance (std::vector< double > &x) |
| void | GetRandomPointMetro (std::vector< double > &x) |
| void | GetRandomPointSamplingMetro (std::vector< double > &x) |
| int | GetNiterationsPerDimension () |
| int | GetNSamplesPer2DBin () |
| int | GetNvar () |
| int | GetNIterationsMax () |
| int | GetNIterations () |
| double | GetRelativePrecision () |
| double | GetError () |
| int | GetNbins () |
| TMinuit * | GetMinuit () |
| int | GetMinuitErrorFlag () |
| TTree * | GetMarkovChainTree () |
| std::vector< double > * | GetMarkovChainPoint () |
| int * | GetMCMCIteration () |
| double * | GetMarkovChainValue () |
| double | GetSAT0 () |
| double | GetSATmin () |
| double * | GetSAP2Temperature () |
| double * | GetSAP2LogProb () |
| std::vector< double > * | GetSAP2x () |
| int * | GetSAP2NIterations () |
Member functions (set) | |
| void | SetMinuitArlist (double *arglist) |
| void | SetFlagIgnorePrevOptimization (bool flag) |
| void | SetParameters (BCParameterSet *par) |
| void | SetVarList (int *varlist) |
| void | SetVar (int index) |
| void | SetIntegrationMethod (BCIntegrate::BCIntegrationMethod method) |
| void | SetMarginalizationMethod (BCIntegrate::BCMarginalizationMethod method) |
| void | SetOptimizationMethod (BCIntegrate::BCOptimizationMethod method) |
| void | SetSASchedule (BCIntegrate::BCSASchedule schedule) |
| void | SetNiterationsPerDimension (int niterations) |
| void | SetNSamplesPer2DBin (int n) |
| void | SetNIterationsMax (int niterations) |
| void | SetRelativePrecision (double relprecision) |
| void | SetNbins (int nbins, int index=-1) |
| void | SetFillErrorBand (bool flag=true) |
| void | UnsetFillErrorBand () |
| void | SetFitFunctionIndexX (int index) |
| void | SetFitFunctionIndexY (int index) |
| void | SetFitFunctionIndices (int indexx, int indexy) |
| void | SetDataPointLowerBoundaries (BCDataPoint *datasetlowerboundaries) |
| void | SetDataPointUpperBoundaries (BCDataPoint *datasetupperboundaries) |
| void | SetDataPointLowerBoundary (int index, double lowerboundary) |
| void | SetDataPointUpperBoundary (int index, double upperboundary) |
| void | WriteMarkovChain (bool flag) |
| void | SetMarkovChainTree (TTree *tree) |
| void | SetMarkovChainInitialPosition (std::vector< double > position) |
| void | SetMarkovChainStepSize (double stepsize) |
| void | SetMarkovChainNIterations (int niterations) |
| void | SetMarkovChainAutoN (bool flag) |
| void | SetMode (std::vector< double > mode) |
| void | SetErrorBandHisto (TH2D *h) |
| void | SetSAT0 (double T0) |
| void | SetSATmin (double Tmin) |
| void | SetFlagWriteSAToFile (bool flag) |
| void | SetSATree (TTree *tree) |
Protected Attributes | |
| int | fNvar |
| int | fNbins |
| int | fNSamplesPer2DBin |
| double | fMarkovChainStepSize |
| int | fMarkovChainNIterations |
| bool | fMarkovChainAutoN |
| BCDataPoint * | fDataPointLowerBoundaries |
| BCDataPoint * | fDataPointUpperBoundaries |
| std::vector< bool > | fDataFixedValues |
| TRandom3 * | fRandom |
| std::vector< double > | fBestFitParameters |
| std::vector< double > | fBestFitParameterErrors |
| std::vector< double > | fBestFitParametersMarginalized |
| std::vector< TH1D * > | fHProb1D |
| std::vector< TH2D * > | fHProb2D |
| bool | fFillErrorBand |
| int | fFitFunctionIndexX |
| int | fFitFunctionIndexY |
| bool | fErrorBandContinuous |
| std::vector< double > | fErrorBandX |
| TH2D * | fErrorBandXY |
| int | fErrorBandNbinsX |
| int | fErrorBandNbinsY |
| TMinuit * | fMinuit |
| double | fMinuitArglist [2] |
| int | fMinuitErrorFlag |
| double | fFlagIgnorePrevOptimization |
| bool | fFlagWriteMarkovChain |
| TTree * | fMarkovChainTree |
| int | fMCMCIteration |
| double | fSAT0 |
| double | fSATmin |
| TTree * | fTreeSA |
| bool | fFlagWriteSAToFile |
| int | fSANIterations |
| double | fSATemperature |
| double | fSALogProb |
| std::vector< double > | fSAx |
Private Member Functions | |
| void | MCMCIterationInterface () |
| void | MCMCFillErrorBand () |
Private Attributes | |
| BCParameterSet * | fx |
| double * | fMin |
| double * | fMax |
| int * | fVarlist |
| int | fNiterPerDimension |
| BCIntegrate::BCIntegrationMethod | fIntegrationMethod |
| BCIntegrate::BCMarginalizationMethod | fMarginalizationMethod |
| BCIntegrate::BCOptimizationMethod | fOptimizationMethod |
| BCIntegrate::BCOptimizationMethod | fOptimizationMethodMode |
| BCIntegrate::BCSASchedule | fSASchedule |
| int | fNIterationsMax |
| int | fNIterations |
| double | fRelativePrecision |
| double | fError |
| int | fNmetro |
| int | fNacceptedMCMC |
| std::vector< double > | fXmetro0 |
| std::vector< double > | fXmetro1 |
| double | fMarkovChainValue |
Definition at line 45 of file BCIntegrate.h.
An enumerator for the integration algorithm
Definition at line 55 of file BCIntegrate.h.
00055 { kIntMonteCarlo, kIntImportance, kIntMetropolis, kIntCuba };
An enumerator for the marginalization algorithm
Definition at line 59 of file BCIntegrate.h.
00059 { kMargMonteCarlo, kMargMetropolis };
An enumerator for the mode finding algorithm
Definition at line 63 of file BCIntegrate.h.
00063 { kOptSA, kOptMetropolis, kOptMinuit };
An enumerator for the Simulated Annealing schedule
Definition at line 67 of file BCIntegrate.h.
00067 { kSACauchy, kSABoltzmann, kSACustom };
| BCIntegrate::BCIntegrate | ( | ) |
The default constructor
Definition at line 31 of file BCIntegrate.cxx.
00031 : BCEngineMCMC() 00032 { 00033 fNvar=0; 00034 fNiterPerDimension = 100; 00035 fNSamplesPer2DBin = 100; 00036 fRandom = new TRandom3(0); 00037 00038 fNIterationsMax = 1000000; 00039 fRelativePrecision = 1e-3; 00040 00041 #ifdef HAVE_CUBA_H 00042 fIntegrationMethod = BCIntegrate::kIntCuba; 00043 #else 00044 fIntegrationMethod = BCIntegrate::kIntMonteCarlo; 00045 #endif 00046 fMarginalizationMethod = BCIntegrate::kMargMetropolis; 00047 fOptimizationMethod = BCIntegrate::kOptMinuit; 00048 fOptimizationMethodMode = BCIntegrate::kOptMinuit; 00049 00050 fNbins=100; 00051 00052 fDataPointLowerBoundaries = 0; 00053 fDataPointUpperBoundaries = 0; 00054 00055 fFillErrorBand = false; 00056 00057 fFitFunctionIndexX = -1; 00058 fFitFunctionIndexY = -1; 00059 00060 fErrorBandNbinsX = 100; 00061 fErrorBandNbinsY = 500; 00062 00063 fErrorBandContinuous = true; 00064 00065 fMinuit = 0; 00066 fMinuitArglist[0] = 20000; 00067 fMinuitArglist[1] = 0.01; 00068 fFlagIgnorePrevOptimization = false; 00069 00070 fFlagWriteMarkovChain = false; 00071 fMarkovChainTree = 0; 00072 00073 fMarkovChainStepSize = 0.1; 00074 00075 fMarkovChainAutoN = true; 00076 00077 fSAT0 = 100.0; 00078 fSATmin = 0.1; 00079 fSASchedule = BCIntegrate::kSACauchy; 00080 fFlagWriteSAToFile = false; 00081 }
| BCIntegrate::BCIntegrate | ( | BCParameterSet * | par | ) |
A constructor
Definition at line 84 of file BCIntegrate.cxx.
00084 : BCEngineMCMC(1) 00085 { 00086 fNvar=0; 00087 fNiterPerDimension = 100; 00088 fNSamplesPer2DBin = 100; 00089 fRandom = new TRandom3(0); 00090 00091 fNbins=100; 00092 00093 this->SetParameters(par); 00094 00095 fDataPointLowerBoundaries = 0; 00096 fDataPointUpperBoundaries = 0; 00097 00098 fFillErrorBand = false; 00099 00100 fFitFunctionIndexX = -1; 00101 fFitFunctionIndexY = -1; 00102 00103 fErrorBandNbinsX = 100; 00104 fErrorBandNbinsY = 200; 00105 00106 fErrorBandContinuous = true; 00107 00108 fMinuit = 0; 00109 fMinuitArglist[0] = 20000; 00110 fMinuitArglist[1] = 0.01; 00111 fFlagIgnorePrevOptimization = false; 00112 00113 fFlagWriteMarkovChain = false; 00114 fMarkovChainTree = 0; 00115 00116 fMarkovChainStepSize = 0.1; 00117 00118 fMarkovChainAutoN = true; 00119 00120 fSAT0 = 100.0; 00121 fSATmin = 0.1; 00122 fTreeSA = false; 00123 }
| BCIntegrate::~BCIntegrate | ( | ) | [virtual] |
The default destructor
Definition at line 126 of file BCIntegrate.cxx.
00127 { 00128 DeleteVarList(); 00129 00130 fx=0; 00131 00132 delete fRandom; 00133 fRandom=0; 00134 00135 if (fMinuit) 00136 delete fMinuit; 00137 00138 int n1 = fHProb1D.size(); 00139 if(n1>0) 00140 { 00141 for (int i=0;i<n1;i++) 00142 delete fHProb1D.at(i); 00143 } 00144 00145 int n2 = fHProb2D.size(); 00146 if(n2>0) 00147 { 00148 for (int i=0;i<n2;i++) 00149 delete fHProb2D.at(i); 00150 } 00151 }
| BCIntegrate::BCIntegrate | ( | ) |
The default constructor
Definition at line 31 of file BCIntegrate.cxx.
00031 : BCEngineMCMC() 00032 { 00033 fNvar=0; 00034 fNiterPerDimension = 100; 00035 fNSamplesPer2DBin = 100; 00036 fRandom = new TRandom3(0); 00037 00038 fNIterationsMax = 1000000; 00039 fRelativePrecision = 1e-3; 00040 00041 #ifdef HAVE_CUBA_H 00042 fIntegrationMethod = BCIntegrate::kIntCuba; 00043 #else 00044 fIntegrationMethod = BCIntegrate::kIntMonteCarlo; 00045 #endif 00046 fMarginalizationMethod = BCIntegrate::kMargMetropolis; 00047 fOptimizationMethod = BCIntegrate::kOptMinuit; 00048 fOptimizationMethodMode = BCIntegrate::kOptMinuit; 00049 00050 fNbins=100; 00051 00052 fDataPointLowerBoundaries = 0; 00053 fDataPointUpperBoundaries = 0; 00054 00055 fFillErrorBand = false; 00056 00057 fFitFunctionIndexX = -1; 00058 fFitFunctionIndexY = -1; 00059 00060 fErrorBandNbinsX = 100; 00061 fErrorBandNbinsY = 500; 00062 00063 fErrorBandContinuous = true; 00064 00065 fMinuit = 0; 00066 fMinuitArglist[0] = 20000; 00067 fMinuitArglist[1] = 0.01; 00068 fFlagIgnorePrevOptimization = false; 00069 00070 fFlagWriteMarkovChain = false; 00071 fMarkovChainTree = 0; 00072 00073 fMarkovChainStepSize = 0.1; 00074 00075 fMarkovChainAutoN = true; 00076 00077 fSAT0 = 100.0; 00078 fSATmin = 0.1; 00079 fSASchedule = BCIntegrate::kSACauchy; 00080 fFlagWriteSAToFile = false; 00081 }
| BCIntegrate::BCIntegrate | ( | BCParameterSet * | par | ) |
A constructor
Definition at line 84 of file BCIntegrate.cxx.
00084 : BCEngineMCMC(1) 00085 { 00086 fNvar=0; 00087 fNiterPerDimension = 100; 00088 fNSamplesPer2DBin = 100; 00089 fRandom = new TRandom3(0); 00090 00091 fNbins=100; 00092 00093 this->SetParameters(par); 00094 00095 fDataPointLowerBoundaries = 0; 00096 fDataPointUpperBoundaries = 0; 00097 00098 fFillErrorBand = false; 00099 00100 fFitFunctionIndexX = -1; 00101 fFitFunctionIndexY = -1; 00102 00103 fErrorBandNbinsX = 100; 00104 fErrorBandNbinsY = 200; 00105 00106 fErrorBandContinuous = true; 00107 00108 fMinuit = 0; 00109 fMinuitArglist[0] = 20000; 00110 fMinuitArglist[1] = 0.01; 00111 fFlagIgnorePrevOptimization = false; 00112 00113 fFlagWriteMarkovChain = false; 00114 fMarkovChainTree = 0; 00115 00116 fMarkovChainStepSize = 0.1; 00117 00118 fMarkovChainAutoN = true; 00119 00120 fSAT0 = 100.0; 00121 fSATmin = 0.1; 00122 fTreeSA = false; 00123 }
| BCIntegrate::~BCIntegrate | ( | ) | [virtual] |
The default destructor
Definition at line 126 of file BCIntegrate.cxx.
00127 { 00128 DeleteVarList(); 00129 00130 fx=0; 00131 00132 delete fRandom; 00133 fRandom=0; 00134 00135 if (fMinuit) 00136 delete fMinuit; 00137 00138 int n1 = fHProb1D.size(); 00139 if(n1>0) 00140 { 00141 for (int i=0;i<n1;i++) 00142 delete fHProb1D.at(i); 00143 } 00144 00145 int n2 = fHProb2D.size(); 00146 if(n2>0) 00147 { 00148 for (int i=0;i<n2;i++) 00149 delete fHProb2D.at(i); 00150 } 00151 }
| BCIntegrate::BCIntegrationMethod BCIntegrate::GetIntegrationMethod | ( | ) | [inline] |
Definition at line 93 of file BCIntegrate.h.
00094 { return fIntegrationMethod; };
| BCIntegrate::BCMarginalizationMethod BCIntegrate::GetMarginalizationMethod | ( | ) | [inline] |
Definition at line 98 of file BCIntegrate.h.
00099 { return fMarginalizationMethod; };
| BCIntegrate::BCOptimizationMethod BCIntegrate::GetOptimizationMethod | ( | ) | [inline] |
Definition at line 103 of file BCIntegrate.h.
00104 { return fOptimizationMethod; };
| BCIntegrate::BCOptimizationMethod BCIntegrate::GetOptimizationMethodMode | ( | ) | [inline] |
Definition at line 108 of file BCIntegrate.h.
00109 { return fOptimizationMethodMode; };
| BCIntegrate::BCSASchedule BCIntegrate::GetSASchedule | ( | ) | [inline] |
Definition at line 113 of file BCIntegrate.h.
00114 { return fSASchedule; };
| void BCIntegrate::GetRandomVector | ( | std::vector< double > & | x | ) |
Fills a vector of random numbers between 0 and 1 into a vector
| A | vector of doubles |
Definition at line 1266 of file BCIntegrate.cxx.
01267 { 01268 double * randx = new double[fNvar]; 01269 01270 fRandom -> RndmArray(fNvar, randx); 01271 01272 for(int i=0;i<fNvar;i++) 01273 x[i] = randx[i]; 01274 01275 delete[] randx; 01276 randx = 0; 01277 }
| void BCIntegrate::GetRandomVectorMetro | ( | std::vector< double > & | x | ) | [virtual] |
Definition at line 1280 of file BCIntegrate.cxx.
01281 { 01282 double * randx = new double[fNvar]; 01283 01284 fRandom -> RndmArray(fNvar, randx); 01285 01286 for(int i=0;i<fNvar;i++) 01287 x[i] = 2.0 * (randx[i] - 0.5) * fMarkovChainStepSize; 01288 01289 delete[] randx; 01290 randx = 0; 01291 }
| double BCIntegrate::GetRandomPoint | ( | std::vector< double > & | x | ) |
Fills a vector of (flat) random numbers in the limits of the parameters and returns the probability at that point
| x | A vector of doubles |
Definition at line 1098 of file BCIntegrate.cxx.
01099 { 01100 this->GetRandomVector(x); 01101 01102 for(int i=0;i<fNvar;i++) 01103 x[i]=fMin[i]+x[i]*(fMax[i]-fMin[i]); 01104 01105 return this->Eval(x); 01106 }
| double BCIntegrate::GetRandomPointImportance | ( | std::vector< double > & | x | ) |
Fills a vector of random numbers in the limits of the parameters sampled by the sampling function and returns the probability at that point
| x | A vector of doubles |
Definition at line 1109 of file BCIntegrate.cxx.
01110 { 01111 double p = 1.1; 01112 double g = 1.0; 01113 01114 while (p>g) 01115 { 01116 this->GetRandomVector(x); 01117 01118 for(int i=0;i<fNvar;i++) 01119 x[i] = fMin[i] + x[i] * (fMax[i]-fMin[i]); 01120 01121 p = fRandom->Rndm(); 01122 01123 g = this -> EvalSampling(x); 01124 } 01125 01126 return this->Eval(x); 01127 }
| void BCIntegrate::GetRandomPointMetro | ( | std::vector< double > & | x | ) |
Fills a vector of random numbers in the limits of the parameters sampled by the probality function and returns the probability at that point (Metropolis)
| x | A vector of doubles |
Definition at line 1154 of file BCIntegrate.cxx.
01155 { 01156 // get new point 01157 this->GetRandomVectorMetro(fXmetro1); 01158 01159 // scale the point to the allowed region and stepsize 01160 int in=1; 01161 for(int i=0;i<fNvar;i++) 01162 { 01163 fXmetro1[i] = fXmetro0[i] + fXmetro1[i] * (fMax[i]-fMin[i]); 01164 01165 // check whether the generated point is inside the allowed region 01166 if( fXmetro1[i]<fMin[i] || fXmetro1[i]>fMax[i] ) 01167 in=0; 01168 } 01169 01170 // calculate the log probabilities and compare old and new point 01171 01172 double p0 = fMarkovChainValue; // old point 01173 double p1 = 0; // new point 01174 int accept=0; 01175 01176 if(in) 01177 { 01178 p1 = this -> LogEval(fXmetro1); 01179 01180 if(p1>=p0) 01181 accept=1; 01182 else 01183 { 01184 double r=log(fRandom->Rndm()); 01185 if(r<p1-p0) 01186 accept=1; 01187 } 01188 } 01189 01190 // fill the return point after the decision 01191 if(accept) 01192 { 01193 // increase counter 01194 fNacceptedMCMC++; 01195 01196 for(int i=0;i<fNvar;i++) 01197 { 01198 fXmetro0[i]=fXmetro1[i]; 01199 x[i]=fXmetro1[i]; 01200 fMarkovChainValue = p1; 01201 } 01202 } 01203 else 01204 for(int i=0;i<fNvar;i++) 01205 { 01206 x[i]=fXmetro0[i]; 01207 fXmetro1[i] = x[i]; 01208 fMarkovChainValue = p0; 01209 } 01210 01211 fNmetro++; 01212 }
| void BCIntegrate::GetRandomPointSamplingMetro | ( | std::vector< double > & | x | ) |
Fills a vector of random numbers in the limits of the parameters sampled by the sampling function and returns the probability at that point (Metropolis)
| x | A vector of doubles |
Definition at line 1215 of file BCIntegrate.cxx.
01216 { 01217 // get new point 01218 this->GetRandomVectorMetro(fXmetro1); 01219 01220 // scale the point to the allowed region and stepsize 01221 int in=1; 01222 for(int i=0;i<fNvar;i++) 01223 { 01224 fXmetro1[i] = fXmetro0[i] + fXmetro1[i] * (fMax[i]-fMin[i]); 01225 01226 // check whether the generated point is inside the allowed region 01227 if( fXmetro1[i]<fMin[i] || fXmetro1[i]>fMax[i] ) 01228 in=0; 01229 } 01230 01231 // calculate the log probabilities and compare old and new point 01232 double p0 = this -> LogEvalSampling(fXmetro0); // old point 01233 double p1=0; // new point 01234 if(in) 01235 p1= this -> LogEvalSampling(fXmetro1); 01236 01237 // compare log probabilities 01238 int accept=0; 01239 if(in) 01240 { 01241 if(p1>=p0) 01242 accept=1; 01243 else 01244 { 01245 double r=log(fRandom->Rndm()); 01246 if(r<p1-p0) 01247 accept=1; 01248 } 01249 } 01250 01251 // fill the return point after the decision 01252 if(accept) 01253 for(int i=0;i<fNvar;i++) 01254 { 01255 fXmetro0[i]=fXmetro1[i]; 01256 x[i]=fXmetro1[i]; 01257 } 01258 else 01259 for(int i=0;i<fNvar;i++) 01260 x[i]=fXmetro0[i]; 01261 01262 fNmetro++; 01263 }
| int BCIntegrate::GetNiterationsPerDimension | ( | ) | [inline] |
Definition at line 151 of file BCIntegrate.h.
00152 { return fNiterPerDimension; };
| int BCIntegrate::GetNSamplesPer2DBin | ( | ) | [inline] |
Definition at line 156 of file BCIntegrate.h.
00157 { return fNSamplesPer2DBin; };
| int BCIntegrate::GetNvar | ( | ) | [inline] |
Definition at line 161 of file BCIntegrate.h.
00162 { return fNvar; };
| int BCIntegrate::GetNIterationsMax | ( | ) | [inline] |
Definition at line 166 of file BCIntegrate.h.
00167 { return fNIterationsMax; };
| int BCIntegrate::GetNIterations | ( | ) | [inline] |
Definition at line 171 of file BCIntegrate.h.
00172 { return fNIterations; };
| double BCIntegrate::GetRelativePrecision | ( | ) | [inline] |
Definition at line 176 of file BCIntegrate.h.
00177 { return fRelativePrecision; };
| double BCIntegrate::GetError | ( | ) | [inline] |
| int BCIntegrate::GetNbins | ( | ) | [inline] |
| TMinuit * BCIntegrate::GetMinuit | ( | ) |
| int BCIntegrate::GetMinuitErrorFlag | ( | ) | [inline] |
| TTree* BCIntegrate::GetMarkovChainTree | ( | ) | [inline] |
| std::vector<double>* BCIntegrate::GetMarkovChainPoint | ( | ) | [inline] |
| int* BCIntegrate::GetMCMCIteration | ( | ) | [inline] |
Returns the iteration of the MCMC
Definition at line 208 of file BCIntegrate.h.
00209 { return &fMCMCIteration; };
| double* BCIntegrate::GetMarkovChainValue | ( | ) | [inline] |
Returns the value of the loglikelihood at the point fXmetro1
Definition at line 213 of file BCIntegrate.h.
00214 { return &fMarkovChainValue; };
| double BCIntegrate::GetSAT0 | ( | ) | [inline] |
Returns the Simulated Annealing starting temperature.
Definition at line 218 of file BCIntegrate.h.
00219 { return fSAT0; }
| double BCIntegrate::GetSATmin | ( | ) | [inline] |
Returns the Simulated Annealing threshhold temperature.
Definition at line 223 of file BCIntegrate.h.
00224 { return fSATmin; }
| double* BCIntegrate::GetSAP2Temperature | ( | ) | [inline] |
| double* BCIntegrate::GetSAP2LogProb | ( | ) | [inline] |
| std::vector<double>* BCIntegrate::GetSAP2x | ( | ) | [inline] |
| int* BCIntegrate::GetSAP2NIterations | ( | ) | [inline] |
| void BCIntegrate::SetMinuitArlist | ( | double * | arglist | ) | [inline] |
Definition at line 243 of file BCIntegrate.h.
00244 { fMinuitArglist[0] = arglist[0]; 00245 fMinuitArglist[1] = arglist[1]; };
| void BCIntegrate::SetFlagIgnorePrevOptimization | ( | bool | flag | ) | [inline] |
| void BCIntegrate::SetParameters | ( | BCParameterSet * | par | ) |
| par | The parameter set which gets translated into array needed for the Monte Carlo integration |
Definition at line 154 of file BCIntegrate.cxx.
00155 { 00156 DeleteVarList(); 00157 00158 fx = par; 00159 fNvar = fx->size(); 00160 fMin = new double[fNvar]; 00161 fMax = new double[fNvar]; 00162 fVarlist = new int[fNvar]; 00163 00164 this->ResetVarlist(1); 00165 00166 for(int i=0;i<fNvar;i++) 00167 { 00168 fMin[i]=(fx->at(i))->GetLowerLimit(); 00169 fMax[i]=(fx->at(i))->GetUpperLimit(); 00170 } 00171 00172 fXmetro0.clear(); 00173 for(int i=0;i<fNvar;i++) 00174 fXmetro0.push_back((fMin[i]+fMax[i])/2.0); 00175 00176 fXmetro1.clear(); 00177 fXmetro1.assign(fNvar,0.); 00178 00179 fMCMCBoundaryMin.clear(); 00180 fMCMCBoundaryMax.clear(); 00181 00182 for(int i=0;i<fNvar;i++) 00183 { 00184 fMCMCBoundaryMin.push_back(fMin[i]); 00185 fMCMCBoundaryMax.push_back(fMax[i]); 00186 fMCMCFlagsFillHistograms.push_back(true); 00187 } 00188 00189 for (int i = int(fMCMCH1NBins.size()); i<fNvar; ++i) 00190 fMCMCH1NBins.push_back(100); 00191 00192 fMCMCNParameters = fNvar; 00193 }
| void BCIntegrate::SetVarList | ( | int * | varlist | ) |
| void BCIntegrate::SetVar | ( | int | index | ) | [inline] |
| index | The index of the variable to be set |
Definition at line 261 of file BCIntegrate.h.
00261 {fVarlist[index]=1;};
| void BCIntegrate::SetIntegrationMethod | ( | BCIntegrate::BCIntegrationMethod | method | ) |
| method | The integration method |
Definition at line 342 of file BCIntegrate.cxx.
00343 { 00344 #ifdef HAVE_CUBA_H 00345 fIntegrationMethod = method; 00346 #else 00347 BCLog::Out(BCLog::warning,BCLog::warning,"!!! This version of BAT is compiled without Cuba."); 00348 BCLog::Out(BCLog::warning,BCLog::warning," Monte Carlo Sampled Mean integration method will be used."); 00349 BCLog::Out(BCLog::warning,BCLog::warning," To be able to use Cuba integration, install Cuba and recompile BAT."); 00350 #endif 00351 }
| void BCIntegrate::SetMarginalizationMethod | ( | BCIntegrate::BCMarginalizationMethod | method | ) | [inline] |
| method | The marginalization method |
Definition at line 269 of file BCIntegrate.h.
00270 { fMarginalizationMethod = method; };
| void BCIntegrate::SetOptimizationMethod | ( | BCIntegrate::BCOptimizationMethod | method | ) | [inline] |
| method | The mode finding method |
Definition at line 274 of file BCIntegrate.h.
00275 { fOptimizationMethod = method; };
| void BCIntegrate::SetSASchedule | ( | BCIntegrate::BCSASchedule | schedule | ) | [inline] |
| method | The Simulated Annealing schedule |
Definition at line 279 of file BCIntegrate.h.
00280 { fSASchedule = schedule; };
| void BCIntegrate::SetNiterationsPerDimension | ( | int | niterations | ) | [inline] |
| niterations | Number of iterations per dimension for Monte Carlo integration. |
Definition at line 284 of file BCIntegrate.h.
00285 { fNiterPerDimension = niterations; };
| void BCIntegrate::SetNSamplesPer2DBin | ( | int | n | ) | [inline] |
| n | Number of samples per 2D bin per variable in the Metropolis marginalization. Default is 100. |
Definition at line 290 of file BCIntegrate.h.
00291 { fNSamplesPer2DBin = n; };
| void BCIntegrate::SetNIterationsMax | ( | int | niterations | ) | [inline] |
| niterations | The maximum number of iterations for Monte Carlo integration |
Definition at line 295 of file BCIntegrate.h.
00296 { fNIterationsMax = niterations; };
| void BCIntegrate::SetRelativePrecision | ( | double | relprecision | ) | [inline] |
| relprecision | The relative precision envisioned for Monte Carlo integration |
Definition at line 301 of file BCIntegrate.h.
00302 { fRelativePrecision = relprecision; };
| void BCIntegrate::SetNbins | ( | int | nbins, | |
| int | index = -1 | |||
| ) |
Set the number of bins for the marginalized distribution of a parameter.
| nbins | Number of bins (default = 100) | |
| index | Index of the parameter. |
Definition at line 226 of file BCIntegrate.cxx.
00227 { 00228 if (fNvar == 0) 00229 return; 00230 00231 // check if index is in range 00232 if (index >= fNvar) 00233 { 00234 BCLog::Out(BCLog::warning, BCLog::warning,"BCIntegrate::SetNbins : Index out of range."); 00235 return; 00236 } 00237 // set for all parameters at once 00238 else if (index < 0) 00239 { 00240 for (int i = 0; i < fNvar; ++i) 00241 this -> SetNbins(nbins, i); 00242 return; 00243 } 00244 00245 // sanity check for nbins 00246 if (nbins <= 0) 00247 nbins = 10; 00248 00249 fMCMCH1NBins[index] = nbins; 00250 00251 return; 00252 00253 // if(n<2) 00254 // { 00255 // BCLog::Out(BCLog::warning, BCLog::warning,"BCIntegrate::SetNbins. Number of bins less than 2 makes no sense. Setting to 2."); 00256 // n=2; 00257 // } 00258 // this -> MCMCSetH1NBins(n, -1); 00259 00260 // fNbins=n; 00261 00262 // fMCMCH1NBins = n; 00263 // fMCMCH2NBinsX = n; 00264 // fMCMCH2NBinsY = n; 00265 }
| void BCIntegrate::SetFillErrorBand | ( | bool | flag = true |
) | [inline] |
| void BCIntegrate::UnsetFillErrorBand | ( | ) | [inline] |
| void BCIntegrate::SetFitFunctionIndexX | ( | int | index | ) | [inline] |
Sets index of the x values in function fits.
| index | Index of the x values |
Definition at line 334 of file BCIntegrate.h.
00335 { fFitFunctionIndexX = index; };
| void BCIntegrate::SetFitFunctionIndexY | ( | int | index | ) | [inline] |
Sets index of the y values in function fits.
| index | Index of the y values |
Definition at line 340 of file BCIntegrate.h.
00341 { fFitFunctionIndexY = index; };
| void BCIntegrate::SetFitFunctionIndices | ( | int | indexx, | |
| int | indexy | |||
| ) | [inline] |
Definition at line 343 of file BCIntegrate.h.
00344 { this -> SetFitFunctionIndexX(indexx); 00345 this -> SetFitFunctionIndexY(indexy); };
| void BCIntegrate::SetDataPointLowerBoundaries | ( | BCDataPoint * | datasetlowerboundaries | ) | [inline] |
Sets the data point containing the lower boundaries of possible data values
Definition at line 350 of file BCIntegrate.h.
00351 { fDataPointLowerBoundaries = datasetlowerboundaries; };
| void BCIntegrate::SetDataPointUpperBoundaries | ( | BCDataPoint * | datasetupperboundaries | ) | [inline] |
Sets the data point containing the upper boundaries of possible data values
Definition at line 356 of file BCIntegrate.h.
00357 { fDataPointUpperBoundaries = datasetupperboundaries; };
| void BCIntegrate::SetDataPointLowerBoundary | ( | int | index, | |
| double | lowerboundary | |||
| ) | [inline] |
Sets the lower boundary of possible data values for a particular variable
Definition at line 362 of file BCIntegrate.h.
00363 { fDataPointLowerBoundaries -> SetValue(index, lowerboundary); };
| void BCIntegrate::SetDataPointUpperBoundary | ( | int | index, | |
| double | upperboundary | |||
| ) | [inline] |
Sets the upper boundary of possible data values for a particular variable
Definition at line 368 of file BCIntegrate.h.
00369 { fDataPointUpperBoundaries -> SetValue(index, upperboundary); };
| void BCIntegrate::WriteMarkovChain | ( | bool | flag | ) | [inline] |
Flag for writing Markov chain to ROOT file (true) or not (false)
Definition at line 374 of file BCIntegrate.h.
00375 { fFlagWriteMarkovChain = flag; fMCMCFlagWriteChainToFile = flag;};
| void BCIntegrate::SetMarkovChainTree | ( | TTree * | tree | ) | [inline] |
Sets the ROOT tree containing the Markov chain
Definition at line 379 of file BCIntegrate.h.
00380 { fMarkovChainTree = tree; };
| void BCIntegrate::SetMarkovChainInitialPosition | ( | std::vector< double > | position | ) |
Sets the initial position for the Markov chain
Definition at line 196 of file BCIntegrate.cxx.
00197 { 00198 int n = position.size(); 00199 00200 fXmetro0.clear(); 00201 00202 for (int i = 0; i < n; ++i) 00203 fXmetro0.push_back(position.at(i)); 00204 }
| void BCIntegrate::SetMarkovChainStepSize | ( | double | stepsize | ) | [inline] |
Sets the step size for Markov chains
Definition at line 388 of file BCIntegrate.h.
00389 { fMarkovChainStepSize = stepsize; };
| void BCIntegrate::SetMarkovChainNIterations | ( | int | niterations | ) | [inline] |
Sets the number of iterations in the markov chain
Definition at line 393 of file BCIntegrate.h.
00394 { fMarkovChainNIterations = niterations; 00395 fMarkovChainAutoN = false; }
| void BCIntegrate::SetMarkovChainAutoN | ( | bool | flag | ) | [inline] |
Sets a flag for automatically calculating the number of iterations
Definition at line 399 of file BCIntegrate.h.
00400 { fMarkovChainAutoN = flag; };
| void BCIntegrate::SetMode | ( | std::vector< double > | mode | ) |
Sets mode
Definition at line 1790 of file BCIntegrate.cxx.
01791 { 01792 if((int)mode.size() == fNvar) 01793 { 01794 fBestFitParameters.clear(); 01795 for (int i = 0; i < fNvar; i++) 01796 fBestFitParameters.push_back(mode[i]); 01797 } 01798 }
| void BCIntegrate::SetErrorBandHisto | ( | TH2D * | h | ) | [inline] |
| void BCIntegrate::SetSAT0 | ( | double | T0 | ) | [inline] |
| T0 | new value for Simulated Annealing starting temperature. |
Definition at line 413 of file BCIntegrate.h.
00414 { fSAT0 = T0; }
| void BCIntegrate::SetSATmin | ( | double | Tmin | ) | [inline] |
| Tmin | new value for Simulated Annealing threshold temperature. |
Definition at line 418 of file BCIntegrate.h.
00419 { fSATmin = Tmin; }
| void BCIntegrate::SetFlagWriteSAToFile | ( | bool | flag | ) | [inline] |
| void BCIntegrate::SetSATree | ( | TTree * | tree | ) | [inline] |
| void BCIntegrate::DeleteVarList | ( | ) |
| void BCIntegrate::ResetVarlist | ( | int | v | ) |
Sets all values of the variable list to a particular value The value
Definition at line 299 of file BCIntegrate.cxx.
| void BCIntegrate::UnsetVar | ( | int | index | ) | [inline] |
Set value of a particular integration variable to 0.
| index | The index of the variable |
Definition at line 446 of file BCIntegrate.h.
00447 { fVarlist[index] = 0; };
| double BCIntegrate::Eval | ( | std::vector< double > | x | ) | [virtual] |
Evaluate the unnormalized probability at a point in parameter space. Method needs to be overloaded by the user.
| x | The point in parameter space |
Reimplemented in BCModel.
Definition at line 306 of file BCIntegrate.cxx.
00307 { 00308 BCLog::Out(BCLog::warning, BCLog::warning, "BCIntegrate::Eval. No function. Function needs to be overloaded."); 00309 return 0; 00310 }
| double BCIntegrate::LogEval | ( | std::vector< double > | x | ) | [virtual] |
Evaluate the natural logarithm of the Eval function. For better numerical stability, this method should also be overloaded by the user.
| x | The point in parameter space |
Reimplemented from BCEngineMCMC.
Reimplemented in BCModel.
Definition at line 313 of file BCIntegrate.cxx.
00314 { 00315 // this method should better also be overloaded 00316 return log(this->Eval(x)); 00317 }
| double BCIntegrate::EvalSampling | ( | std::vector< double > | x | ) | [virtual] |
Evaluate the sampling function at a point in parameter space. Method needs to be overloaded by the user.
| x | The point in parameter space |
Reimplemented in BCModel.
Definition at line 320 of file BCIntegrate.cxx.
00321 { 00322 BCLog::Out(BCLog::warning, BCLog::warning, "BCIntegrate::EvalSampling. No function. Function needs to be overloaded."); 00323 return 0; 00324 }
| double BCIntegrate::LogEvalSampling | ( | std::vector< double > | x | ) |
Evaluate the natural logarithm of the EvalSampling function. Method needs to be overloaded by the user.
| x | The point in parameter space |
Definition at line 327 of file BCIntegrate.cxx.
| double BCIntegrate::EvalPrint | ( | std::vector< double > | x | ) |
Evaluate the un-normalized probability at a point in parameter space and prints the result to the log.
| x | The point in parameter space |
Definition at line 333 of file BCIntegrate.cxx.
00334 { 00335 double val=this->Eval(x); 00336 BCLog::Out(BCLog::detail, BCLog::detail, Form("BCIntegrate::EvalPrint. Value: %d.", val)); 00337 00338 return val; 00339 }
| virtual double BCIntegrate::FitFunction | ( | std::vector< double > | x, | |
| std::vector< double > | parameters | |||
| ) | [inline, virtual] |
Defines a fit function.
| parameters | A set of parameter values | |
| x | A vector of x-values |
Reimplemented in BCEfficiencyFitter, BCGraphFitter, and BCHistogramFitter.
Definition at line 490 of file BCIntegrate.h.
| double BCIntegrate::Integrate | ( | ) |
Does the integration over the un-normalized probability.
Definition at line 354 of file BCIntegrate.cxx.
00355 { 00356 std::vector <double> parameter; 00357 parameter.assign(fNvar, 0.0); 00358 00359 switch(fIntegrationMethod) 00360 { 00361 case BCIntegrate::kIntMonteCarlo: 00362 return IntegralMC(parameter); 00363 00364 case BCIntegrate::kIntMetropolis: 00365 return this -> IntegralMetro(parameter); 00366 00367 case BCIntegrate::kIntImportance: 00368 return this -> IntegralImportance(parameter); 00369 00370 case BCIntegrate::kIntCuba: 00371 #ifdef HAVE_CUBA_H 00372 return this -> CubaIntegrate(); 00373 #else 00374 BCLog::Out(BCLog::error,BCLog::error,"!!! This version of BAT is compiled without Cuba."); 00375 BCLog::Out(BCLog::error,BCLog::error," Use other integration methods or install Cuba and recompile BAT."); 00376 #endif 00377 } 00378 00379 BCLog::Out(BCLog::error, BCLog::error, 00380 Form("BCIntegrate::Integrate : Invalid integration method: %d", fIntegrationMethod)); 00381 00382 return 0; 00383 }
| double BCIntegrate::IntegralMC | ( | std::vector< double > | x, | |
| int * | varlist | |||
| ) |
Perfoms the Monte Carlo integration. For details see documentation.
| x | An initial point in parameter space | |
| varlist | A list of variables |
Definition at line 386 of file BCIntegrate.cxx.
00387 { 00388 this->SetVarList(varlist); 00389 return IntegralMC(x); 00390 }
| double BCIntegrate::IntegralMC | ( | std::vector< double > | x | ) |
Definition at line 393 of file BCIntegrate.cxx.
00394 { 00395 // count the variables to integrate over 00396 int NvarNow=0; 00397 00398 for(int i = 0; i < fNvar; i++) 00399 if(fVarlist[i]) 00400 NvarNow++; 00401 00402 // print to log 00403 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Running MC integation over %i dimensions.", NvarNow)); 00404 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Maximum number of iterations: %i", this->GetNIterationsMax())); 00405 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Aimed relative precision: %e", this->GetRelativePrecision())); 00406 00407 // calculate D (the integrated volume) 00408 double D = 1.0; 00409 for(int i = 0; i < fNvar; i++) 00410 if (fVarlist[i]) 00411 D = D * (fMax[i] - fMin[i]); 00412 00413 // reset variables 00414 double pmax = 0.0; 00415 double sumW = 0.0; 00416 double sumW2 = 0.0; 00417 double integral = 0.0; 00418 double variance = 0.0; 00419 double relprecision = 1.0; 00420 00421 std::vector <double> randx; 00422 randx.assign(fNvar, 0.0); 00423 00424 // reset number of iterations 00425 fNIterations = 0; 00426 00427 // iterate while precision is not reached and the number of iterations is lower than maximum number of iterations 00428 while ((fRelativePrecision < relprecision && fNIterationsMax > fNIterations) || fNIterations < 10) 00429 { 00430 // increase number of iterations 00431 fNIterations++; 00432 00433 // get random numbers 00434 this -> GetRandomVector(randx); 00435 00436 // scale random numbers 00437 for(int i = 0; i < fNvar; i++) 00438 { 00439 if(fVarlist[i]) 00440 randx[i]=fMin[i]+randx[i]*(fMax[i]-fMin[i]); 00441 else 00442 randx[i]=x[i]; 00443 } 00444 00445 // evaluate function at sampled point 00446 double value = this->Eval(randx); 00447 00448 // add value to sum and sum of squares 00449 sumW += value; 00450 sumW2 += value * value; 00451 00452 // search for maximum probability 00453 if (value > pmax) 00454 { 00455 // set new maximum value 00456 pmax = value; 00457 00458 // delete old best fit parameter values 00459 fBestFitParameters.clear(); 00460 00461 // write best fit parameters 00462 for(int i = 0; i < fNvar; i++) 00463 fBestFitParameters.push_back(randx.at(i)); 00464 } 00465 00466 // calculate integral and variance 00467 integral = D * sumW / fNIterations; 00468 00469 if (fNIterations%10000 == 0) 00470 { 00471 variance = (1.0 / double(fNIterations)) * (D * D * sumW2 / double(fNIterations) - integral * integral); 00472 double error = sqrt(variance); 00473 relprecision = error / integral; 00474 00475 BCLog::Out(BCLog::debug, BCLog::debug, 00476 Form("BCIntegrate::IntegralMC. Iteration %i, integral: %e +- %e.", fNIterations, integral, error)); 00477 } 00478 } 00479 00480 fError = variance / fNIterations; 00481 00482 // print to log 00483 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Result of integration: %e +- %e in %i iterations.", integral, sqrt(variance), fNIterations)); 00484 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Obtained relative precision: %e. ", sqrt(variance) / integral)); 00485 00486 return integral; 00487 }
| double BCIntegrate::IntegralMetro | ( | std::vector< double > | x | ) |
Perfoms the Metropolis Monte Carlo integration. For details see documentation.
| x | An initial point in parameter space |
Definition at line 491 of file BCIntegrate.cxx.
00492 { 00493 // print debug information 00494 BCLog::Out(BCLog::debug, BCLog::debug, Form("BCIntegrate::IntegralMetro. Integate over %i dimensions.", fNvar)); 00495 00496 // get total number of iterations 00497 double Niter = pow(fNiterPerDimension, fNvar); 00498 00499 // print if more than 100,000 iterations 00500 if(Niter>1e5) 00501 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Total number of iterations in Metropolis: %d.", Niter)); 00502 00503 // reset sum 00504 double sumI = 0; 00505 00506 // prepare Metropolis algorithm 00507 std::vector <double> randx; 00508 randx.assign(fNvar,0.); 00509 InitMetro(); 00510 00511 // prepare maximum value 00512 double pmax = 0.0; 00513 00514 // loop over iterations 00515 for(int i = 0; i <= Niter; i++) 00516 { 00517 // get random point from sampling distribution 00518 this -> GetRandomPointSamplingMetro(randx); 00519 00520 // calculate probability at random point 00521 double val_f = this -> Eval(randx); 00522 00523 // calculate sampling distributions at that point 00524 double val_g = this -> EvalSampling(randx); 00525 00526 // add ratio to sum 00527 if (val_g > 0) 00528 sumI += val_f / val_g; 00529 00530 // search for maximum probability 00531 if (val_f > pmax) 00532 { 00533 // set new maximum value 00534 pmax = val_f; 00535 00536 // delete old best fit parameter values 00537 fBestFitParameters.clear(); 00538 00539 // write best fit parameters 00540 for(int i = 0; i < fNvar; i++) 00541 fBestFitParameters.push_back(randx.at(i)); 00542 } 00543 00544 // write intermediate results 00545 if((i+1)%100000 == 0) 00546 BCLog::Out(BCLog::debug, BCLog::debug, 00547 Form("BCIntegrate::IntegralMetro. Iteration %d, integral: %d.", i+1, sumI/(i+1))); 00548 } 00549 00550 // calculate integral 00551 double result=sumI/Niter; 00552 00553 // print debug information 00554 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Integral is %f in %i iterations.", result, Niter)); 00555 00556 return result; 00557 }
| double BCIntegrate::IntegralImportance | ( | std::vector< double > | x | ) |
Perfoms the importance sampling Monte Carlo integration. For details see documentation.
| x | An initial point in parameter space |
Definition at line 560 of file BCIntegrate.cxx.
00561 { 00562 // print debug information 00563 BCLog::Out(BCLog::debug, BCLog::debug, Form("BCIntegrate::IntegralImportance. Integate over %i dimensions.", fNvar)); 00564 00565 // get total number of iterations 00566 double Niter = pow(fNiterPerDimension, fNvar); 00567 00568 // print if more than 100,000 iterations 00569 if(Niter>1e5) 00570 BCLog::Out(BCLog::detail, BCLog::detail, Form("BCIntegrate::IntegralImportance. Total number of iterations: %d.", Niter)); 00571 00572 // reset sum 00573 double sumI = 0; 00574 00575 std::vector <double> randx; 00576 randx.assign(fNvar,0.); 00577 00578 // prepare maximum value 00579 double pmax = 0.0; 00580 00581 // loop over iterations 00582 for(int i = 0; i <= Niter; i++) 00583 { 00584 // get random point from sampling distribution 00585 this -> GetRandomPointImportance(randx); 00586 00587 // calculate probability at random point 00588 double val_f = this -> Eval(randx); 00589 00590 // calculate sampling distributions at that point 00591 double val_g = this -> EvalSampling(randx); 00592 00593 // add ratio to sum 00594 if (val_g > 0) 00595 sumI += val_f / val_g; 00596 00597 // search for maximum probability 00598 if (val_f > pmax) 00599 { 00600 // set new maximum value 00601 pmax = val_f; 00602 00603 // delete old best fit parameter values 00604 fBestFitParameters.clear(); 00605 00606 // write best fit parameters 00607 for(int i = 0; i < fNvar; i++) 00608 fBestFitParameters.push_back(randx.at(i)); 00609 } 00610 00611 // write intermediate results 00612 if((i+1)%100000 == 0) 00613 BCLog::Out(BCLog::debug, BCLog::debug, 00614 Form("BCIntegrate::IntegralImportance. Iteration %d, integral: %d.", i+1, sumI/(i+1))); 00615 } 00616 00617 // calculate integral 00618 double result=sumI/Niter; 00619 00620 // print debug information 00621 BCLog::Out(BCLog::debug, BCLog::debug, Form("BCIntegrate::IntegralImportance. Integral %f in %i iterations.", result, Niter)); 00622 00623 return result; 00624 }
| double BCIntegrate::CubaIntegrate | ( | int | method, | |
| std::vector< double > | parameters_double, | |||
| std::vector< double > | parameters_int | |||
| ) |
Calculate integral using the Cuba library. For details see documentation.
| method | A short cut for the method | |
| parameters_double | A vector of parameters (double) | |
| parameters_int | A vector of parameters (int) |
Definition at line 1938 of file BCIntegrate.cxx.
01939 { 01940 #ifdef HAVE_CUBA_H 01941 const int NDIM = int(fx ->size()); 01942 const int NCOMP = 1; 01943 01944 const double EPSREL = parameters_double[0]; 01945 const double EPSABS = parameters_double[1]; 01946 const int VERBOSE = int(parameters_int[0]); 01947 const int MINEVAL = int(parameters_int[1]); 01948 const int MAXEVAL = int(parameters_int[2]); 01949 01950 int neval; 01951 int fail; 01952 int nregions; 01953 double integral[NCOMP]; 01954 double error[NCOMP]; 01955 double prob[NCOMP]; 01956 01957 global_this = this; 01958 01959 integrand_t an_integrand = &BCIntegrate::CubaIntegrand; 01960 01961 if (method == 0) 01962 { 01963 // set VEGAS specific parameters 01964 const int NSTART = int(parameters_int[3]); 01965 const int NINCREASE = int(parameters_int[4]); 01966 01967 // call VEGAS integration method 01968 Vegas(NDIM, NCOMP, an_integrand, 01969 EPSREL, EPSABS, VERBOSE, MINEVAL, MAXEVAL, 01970 NSTART, NINCREASE, 01971 &neval, &fail, integral, error, prob); 01972 } 01973 else if (method == 1) 01974 { 01975 // set SUAVE specific parameters 01976 const int LAST = int(parameters_int[3]); 01977 const int NNEW = int(parameters_int[4]); 01978 const int FLATNESS = int(parameters_int[5]); 01979 01980 // call SUAVE integration method 01981 Suave(NDIM, NCOMP, an_integrand, 01982 EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, 01983 NNEW, FLATNESS, 01984 &nregions, &neval, &fail, integral, error, prob); 01985 } 01986 else if (method == 2) 01987 { 01988 // set DIVONNE specific parameters 01989 const int KEY1 = int(parameters_int[3]); 01990 const int KEY2 = int(parameters_int[4]); 01991 const int KEY3 = int(parameters_int[5]); 01992 const int MAXPASS = int(parameters_int[6]); 01993 const int BORDER = int(parameters_int[7]); 01994 const int MAXCHISQ = int(parameters_int[8]); 01995 const int MINDEVIATION = int(parameters_int[9]); 01996 const int NGIVEN = int(parameters_int[10]); 01997 const int LDXGIVEN = int(parameters_int[11]); 01998 const int NEXTRA = int(parameters_int[12]); 01999 02000 // call DIVONNE integration method 02001 Divonne(NDIM, NCOMP, an_integrand, 02002 EPSREL, EPSABS, VERBOSE, MINEVAL, MAXEVAL, 02003 KEY1, KEY2, KEY3, MAXPASS, BORDER, MAXCHISQ, MINDEVIATION, 02004 NGIVEN, LDXGIVEN, NULL, NEXTRA, NULL, 02005 &nregions, &neval, &fail, integral, error, prob); 02006 } 02007 else if (method == 3) 02008 { 02009 // set CUHRE specific parameters 02010 const int LAST = int(parameters_int[3]); 02011 const int KEY = int(parameters_int[4]); 02012 02013 // call CUHRE integration method 02014 Cuhre(NDIM, NCOMP, an_integrand, 02015 EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, KEY, 02016 &nregions, &neval, &fail, integral, error, prob); 02017 } 02018 else 02019 { 02020 std::cout << " Integration method not available. " << std::endl; 02021 integral[0] = -1e99; 02022 } 02023 02024 if (fail != 0) 02025 { 02026 std::cout << " Warning, integral did not converge with the given set of parameters. "<< std::endl; 02027 std::cout << " neval = " << neval << std::endl; 02028 std::cout << " fail = " << fail << std::endl; 02029 std::cout << " integral = " << integral[0] << std::endl; 02030 std::cout << " error = " << error[0] << std::endl; 02031 std::cout << " prob = " << prob[0] << std::endl; 02032 } 02033 02034 return integral[0] / 1e99; 02035 #else 02036 BCLog::Out(BCLog::error,BCLog::error,"!!! This version of BAT is compiled without Cuba."); 02037 BCLog::Out(BCLog::error,BCLog::error," Use other integration methods or install Cuba and recompile BAT."); 02038 return 0.; 02039 #endif 02040 }
| double BCIntegrate::CubaIntegrate | ( | ) |
Definition at line 1901 of file BCIntegrate.cxx.
01902 { 01903 #ifdef HAVE_CUBA_H 01904 double EPSREL = 1e-3; 01905 double EPSABS = 1e-12; 01906 double VERBOSE = 0; 01907 double MINEVAL = 0; 01908 double MAXEVAL = 2000000; 01909 double NSTART = 25000; 01910 double NINCREASE = 25000; 01911 01912 std::vector<double> parameters_double; 01913 std::vector<double> parameters_int; 01914 01915 parameters_double.push_back(EPSREL); 01916 parameters_double.push_back(EPSABS); 01917 01918 parameters_int.push_back(VERBOSE); 01919 parameters_int.push_back(MINEVAL); 01920 parameters_int.push_back(MAXEVAL); 01921 parameters_int.push_back(NSTART); 01922 parameters_int.push_back(NINCREASE); 01923 01924 // print to log 01925 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Running Cuba/Vegas integation over %i dimensions.", fNvar)); 01926 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Maximum number of iterations: %i", (int)MAXEVAL)); 01927 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Aimed relative precision: %e", EPSREL)); 01928 01929 return this -> CubaIntegrate(0, parameters_double, parameters_int); 01930 #else 01931 BCLog::Out(BCLog::error,BCLog::error,"!!! This version of BAT is compiled without Cuba."); 01932 BCLog::Out(BCLog::error,BCLog::error," Use other integration methods or install Cuba and recompile BAT."); 01933 return 0.; 01934 #endif 01935 }
| void BCIntegrate::CubaIntegrand | ( | const int * | ndim, | |
| const double | xx[], | |||
| const int * | ncomp, | |||
| double | ff[] | |||
| ) | [static] |
Integrand for the Cuba library.
| ndim | The number of dimensions to integrate over | |
| xx | The point in parameter space to integrate over (scaled to 0 - 1 per dimension) | |
| ncomp | The number of components of the integrand (usually 1) | |
| ff | The function value |
Definition at line 1863 of file BCIntegrate.cxx.
01865 { 01866 #ifdef HAVE_CUBA_H 01867 // scale variables 01868 double jacobian = 1.0; 01869 01870 std::vector<double> scaled_parameters; 01871 01872 for (int i = 0; i < *ndim; i++) 01873 { 01874 double range = global_this -> fx -> at(i) -> GetUpperLimit() - global_this -> fx -> at(i) -> GetLowerLimit(); 01875 01876 // multiply range to jacobian 01877 jacobian *= range; 01878 01879 // get the scaled parameter value 01880 scaled_parameters.push_back(global_this -> fx -> at(i) -> GetLowerLimit() + xx[i] * range); 01881 } 01882 01883 // call function to integrate 01884 ff[0] = global_this -> Eval(scaled_parameters); 01885 01886 // multiply jacobian 01887 ff[0] *= jacobian; 01888 01889 // multiply fudge factor 01890 ff[0] *= 1e99; 01891 01892 // remove parameter vector 01893 scaled_parameters.clear(); 01894 #else 01895 BCLog::Out(BCLog::error,BCLog::error,"!!! This version of BAT is compiled without Cuba."); 01896 BCLog::Out(BCLog::error,BCLog::error," Use other integration methods or install Cuba and recompile BAT."); 01897 #endif 01898 }
| TH1D * BCIntegrate::Marginalize | ( | BCParameter * | parameter | ) |
Performs the marginalization with respect to one parameter.
| parameter | The parameter w.r.t. which the marginalization is performed |
Definition at line 627 of file BCIntegrate.cxx.
00628 { 00629 BCLog::Out(BCLog::detail, BCLog::detail, 00630 Form(" --> Marginalizing model wrt. parameter %s using method %d.", parameter->GetName().data(), fMarginalizationMethod)); 00631 00632 switch(fMarginalizationMethod) 00633 { 00634 case BCIntegrate::kMargMonteCarlo: 00635 return MarginalizeByIntegrate(parameter); 00636 00637 case BCIntegrate::kMargMetropolis: 00638 return MarginalizeByMetro(parameter); 00639 } 00640 00641 BCLog::Out(BCLog::warning, BCLog::warning, 00642 Form("BCIntegrate::Marginalize. Invalid marginalization method: %d. Return 0.", fMarginalizationMethod)); 00643 00644 return 0; 00645 }
| TH2D * BCIntegrate::Marginalize | ( | BCParameter * | parameter1, | |
| BCParameter * | parameter2 | |||
| ) |
Performs the marginalization with respect to two parameters.
| parameter1 | The first parameter w.r.t. which the marginalization is performed | |
| parameter2 | The second parameter w.r.t. which the marginalization is performed |
Definition at line 648 of file BCIntegrate.cxx.
00649 { 00650 switch(fMarginalizationMethod) 00651 { 00652 case BCIntegrate::kMargMonteCarlo: 00653 return MarginalizeByIntegrate(parameter1,parameter2); 00654 00655 case BCIntegrate::kMargMetropolis: 00656 return MarginalizeByMetro(parameter1,parameter2); 00657 } 00658 00659 BCLog::Out(BCLog::warning, BCLog::warning, 00660 Form("BCIntegrate::Marginalize. Invalid marginalization method: %d. Return 0.", fMarginalizationMethod)); 00661 00662 return 0; 00663 }
| TH1D * BCIntegrate::MarginalizeByIntegrate | ( | BCParameter * | parameter | ) |
Performs the marginalization with respect to one parameter using the simple Monte Carlo technique.
| parameter | The parameter w.r.t. which the marginalization is performed |
Definition at line 666 of file BCIntegrate.cxx.
00667 { 00668 // set parameter to marginalize 00669 this->ResetVarlist(1); 00670 int index = parameter->GetIndex(); 00671 this->UnsetVar(index); 00672 00673 // define histogram 00674 double hmin = parameter -> GetLowerLimit(); 00675 double hmax = parameter -> GetUpperLimit(); 00676 double hdx = (hmax - hmin) / double(fNbins); 00677 TH1D * hist = new TH1D("hist","", fNbins, hmin, hmax); 00678 00679 // integrate 00680 std::vector <double> randx; 00681 randx.assign(fNvar, 0.0); 00682 00683 for(int i=0;i<=fNbins;i++) 00684 { 00685 randx[index] = hmin + (double)i * hdx; 00686 00687 double val = IntegralMC(randx); 00688 hist->Fill(randx[index], val); 00689 } 00690 00691 // normalize 00692 hist -> Scale( 1./hist->Integral("width") ); 00693 00694 return hist; 00695 }
| TH2D * BCIntegrate::MarginalizeByIntegrate | ( | BCParameter * | parameter1, | |
| BCParameter * | parameter2 | |||
| ) |
Performs the marginalization with respect to two parameters using the simple Monte Carlo technique.
| parameter1 | The first parameter w.r.t. which the marginalization is performed | |
| parameter2 | The second parameter w.r.t. which the marginalization is performed |
Definition at line 698 of file BCIntegrate.cxx.
00699 { 00700 // set parameter to marginalize 00701 this->ResetVarlist(1); 00702 int index1 = parameter1->GetIndex(); 00703 this->UnsetVar(index1); 00704 int index2 = parameter2->GetIndex(); 00705 this->UnsetVar(index2); 00706 00707 // define histogram 00708 double hmin1 = parameter1 -> GetLowerLimit(); 00709 double hmax1 = parameter1 -> GetUpperLimit(); 00710 double hdx1 = (hmax1 - hmin1) / double(fNbins); 00711 00712 double hmin2 = parameter2 -> GetLowerLimit(); 00713 double hmax2 = parameter2 -> GetUpperLimit(); 00714 double hdx2 = (hmax2 - hmin2) / double(fNbins); 00715 00716 TH2D * hist = new TH2D(Form("hist_%s_%s", parameter1 -> GetName().data(), parameter2 -> GetName().data()),"", 00717 fNbins, hmin1, hmax1, 00718 fNbins, hmin2, hmax2); 00719 00720 // integrate 00721 std::vector <double> randx; 00722 randx.assign(fNvar, 0.0); 00723 00724 for(int i=0;i<=fNbins;i++) 00725 { 00726 randx[index1] = hmin1 + (double)i * hdx1; 00727 for(int j=0;j<=fNbins;j++) 00728 { 00729 randx[index2] = hmin2 + (double)j * hdx2; 00730 00731 double val = IntegralMC(randx); 00732 hist->Fill(randx[index1],randx[index2], val); 00733 } 00734 } 00735 00736 // normalize 00737 hist -> Scale(1.0/hist->Integral("width")); 00738 00739 return hist; 00740 }
| TH1D * BCIntegrate::MarginalizeByMetro | ( | BCParameter * | parameter | ) |
Performs the marginalization with respect to one parameter using the Metropolis algorithm.
| parameter | The parameter w.r.t. which the marginalization is performed |
Definition at line 743 of file BCIntegrate.cxx.
00744 { 00745 int niter = fMarkovChainNIterations; 00746 00747 if (fMarkovChainAutoN == true) 00748 niter = fNbins*fNbins*fNSamplesPer2DBin*fNvar; 00749 00750 BCLog::Out(BCLog::detail, BCLog::detail, 00751 Form(" --> Number of samples in Metropolis marginalization: %d.", niter)); 00752 00753 // set parameter to marginalize 00754 int index = parameter->GetIndex(); 00755 00756 // define histogram 00757 double hmin = parameter -> GetLowerLimit(); 00758 double hmax = parameter -> GetUpperLimit(); 00759 TH1D * hist = new TH1D("hist","", fNbins, hmin, hmax); 00760 00761 // prepare Metro 00762 std::vector <double> randx; 00763 randx.assign(fNvar, 0.0); 00764 InitMetro(); 00765 00766 for(int i=0;i<=niter;i++) 00767 { 00768 GetRandomPointMetro(randx); 00769 hist->Fill(randx[index]); 00770 } 00771 00772 // normalize 00773 hist -> Scale(1.0/hist->Integral("width")); 00774 00775 return hist; 00776 }
| TH2D * BCIntegrate::MarginalizeByMetro | ( | BCParameter * | parameter1, | |
| BCParameter * | parameter2 | |||
| ) |
Performs the marginalization with respect to two parameters using the Metropolis algorithm.
| parameter1 | The first parameter w.r.t. which the marginalization is performed | |
| parameter2 | The second parameter w.r.t. which the marginalization is performed |
Definition at line 779 of file BCIntegrate.cxx.
00780 { 00781 int niter=fNbins*fNbins*fNSamplesPer2DBin*fNvar; 00782 00783 // set parameter to marginalize 00784 int index1 = parameter1->GetIndex(); 00785 int index2 = parameter2->GetIndex(); 00786 00787 // define histogram 00788 double hmin1 = parameter1 -> GetLowerLimit(); 00789 double hmax1 = parameter1 -> GetUpperLimit(); 00790 00791 double hmin2 = parameter2 -> GetLowerLimit(); 00792 double hmax2 = parameter2 -> GetUpperLimit(); 00793 00794 TH2D * hist = new TH2D(Form("hist_%s_%s", parameter1 -> GetName().data(), parameter2 -> GetName().data()),"", 00795 fNbins, hmin1, hmax1, 00796 fNbins, hmin2, hmax2); 00797 00798 // prepare Metro 00799 std::vector <double> randx; 00800 randx.assign(fNvar, 0.0); 00801 InitMetro(); 00802 00803 for(int i=0;i<=niter;i++) 00804 { 00805 GetRandomPointMetro(randx); 00806 hist->Fill(randx[index1],randx[index2]); 00807 } 00808 00809 // normalize 00810 hist -> Scale(1.0/hist->Integral("width")); 00811 00812 return hist; 00813 }
| int BCIntegrate::MarginalizeAllByMetro | ( | const char * | name | ) |
Performs the marginalization with respect to every single parameter as well as with respect all combinations to two parameters using the Metropolis algorithm.
| name | Basename for the histograms (e.g. model name) |
Definition at line 816 of file BCIntegrate.cxx.
00817 { 00818 int niter=fNbins*fNbins*fNSamplesPer2DBin*fNvar; 00819 00820 BCLog::Out(BCLog::detail, BCLog::detail, 00821 Form(" --> Number of samples in Metropolis marginalization: %d.", niter)); 00822 00823 // define 1D histograms 00824 for(int i=0;i<fNvar;i++) 00825 { 00826 double hmin1 = fx->at(i) -> GetLowerLimit(); 00827 double hmax1 = fx->at(i) -> GetUpperLimit(); 00828 00829 TH1D * h1 = new TH1D( 00830 TString::Format("h%s_%s", name, fx->at(i) -> GetName().data()),"", 00831 fNbins, hmin1, hmax1); 00832 00833 fHProb1D.push_back(h1); 00834 } 00835 00836 // define 2D histograms 00837 for(int i=0;i<fNvar-1;i++) 00838 for(int j=i+1;j<fNvar;j++) 00839 { 00840 double hmin1 = fx->at(i) -> GetLowerLimit(); 00841 double hmax1 = fx->at(i) -> GetUpperLimit(); 00842 00843 double hmin2 = fx->at(j) -> GetLowerLimit(); 00844 double hmax2 = fx->at(j) -> GetUpperLimit(); 00845 00846 TH2D * h2 = new TH2D( 00847 TString::Format("h%s_%s_%s", name, fx->at(i) -> GetName().data(), fx->at(j) -> GetName().data()),"", 00848 fNbins, hmin1, hmax1, 00849 fNbins, hmin2, hmax2); 00850 00851 fHProb2D.push_back(h2); 00852 } 00853 00854 // get number of 2d distributions 00855 int nh2d = fHProb2D.size(); 00856 00857 BCLog::Out(BCLog::detail, BCLog::detail, 00858 Form(" --> Marginalizing %d 1D distributions and %d 2D distributions.", fNvar, nh2d)); 00859 00860 // prepare function fitting 00861 double dx = 0.0; 00862 double dy = 0.0; 00863 00864 if (fFitFunctionIndexX >= 0) 00865 { 00866 dx = (fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexX) - 00867 fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexX)) 00868 / double(fErrorBandNbinsX); 00869 00870 dx = (fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexY) - 00871 fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexY)) 00872 / double(fErrorBandNbinsY); 00873 00874 fErrorBandXY = new TH2D( 00875 TString::Format("errorbandxy_%d",BCLog::GetHIndex()), "", 00876 fErrorBandNbinsX, 00877 fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexX) - 0.5 * dx, 00878 fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexX) + 0.5 * dx, 00879 fErrorBandNbinsY, 00880 fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexY) - 0.5 * dy, 00881 fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexY) + 0.5 * dy); 00882 fErrorBandXY -> SetStats(kFALSE); 00883 00884 for (int ix = 1; ix <= fErrorBandNbinsX; ++ix) 00885 for (int iy = 1; iy <= fErrorBandNbinsX; ++iy) 00886 fErrorBandXY -> SetBinContent(ix, iy, 0.0); 00887 } 00888 00889 // prepare Metro 00890 std::vector <double> randx; 00891 00892 randx.assign(fNvar, 0.0); 00893 InitMetro(); 00894 00895 // reset counter 00896 fNacceptedMCMC = 0; 00897 00898 // run Metro and fill histograms 00899 for(int i=0;i<=niter;i++) 00900 { 00901 GetRandomPointMetro(randx); 00902 00903 // save this point to the markov chain in the ROOT file 00904 if (fFlagWriteMarkovChain) 00905 { 00906 fMCMCIteration = i; 00907 fMarkovChainTree -> Fill(); 00908 } 00909 00910 for(int j=0;j<fNvar;j++) 00911 fHProb1D[j] -> Fill( randx[j] ); 00912 00913 int ih=0; 00914 for(int j=0;j<fNvar-1;j++) 00915 for(int k=j+1;k<fNvar;k++) 00916 { 00917 fHProb2D[ih] -> Fill(randx[j],randx[k]); 00918 ih++; 00919 } 00920 00921 if((i+1)%100000==0) 00922 BCLog::Out(BCLog::debug, BCLog::debug, 00923 Form("BCIntegrate::MarginalizeAllByMetro. %d samples done.",i+1)); 00924 00925 // function fitting 00926 00927 if (fFitFunctionIndexX >= 0) 00928 { 00929 // loop over all possible x values ... 00930 if (fErrorBandContinuous) 00931 { 00932 double x = 0; 00933 00934 for (int ix = 0; ix < fErrorBandNbinsX; ix++) 00935 { 00936 // calculate x 00937 x = fErrorBandXY -> GetXaxis() -> GetBinCenter(ix + 1); 00938 00939 // calculate y 00940 std::vector <double> xvec; 00941 xvec.push_back(x); 00942 double y = this -> FitFunction(xvec, randx); 00943 xvec.clear(); 00944 00945 // fill histogram 00946 fErrorBandXY -> Fill(x, y); 00947 } 00948 } 00949 00950 // ... or evaluate at the data point x-values 00951 else 00952 { 00953 int ndatapoints = int(fErrorBandX.size()); 00954 double x = 0; 00955 00956 for (int ix = 0; ix < ndatapoints; ++ix) 00957 { 00958 // calculate x 00959 x = fErrorBandX.at(ix); 00960 00961 // calculate y 00962 std::vector <double> xvec; 00963 xvec.push_back(x); 00964 double y = this -> FitFunction(xvec, randx); 00965 xvec.clear(); 00966 00967 // fill histogram 00968 fErrorBandXY -> Fill(x, y); 00969 } 00970 } 00971 } 00972 } 00973 00974 // normalize histograms 00975 for(int i=0;i<fNvar;i++) 00976 fHProb1D[i] -> Scale( 1./fHProb1D[i]->Integral("width") ); 00977 00978 for (int i=0;i<nh2d;i++) 00979 fHProb2D[i] -> Scale( 1./fHProb2D[i]->Integral("width") ); 00980 00981 if (fFitFunctionIndexX >= 0) 00982 fErrorBandXY -> Scale(1.0/fErrorBandXY -> Integral() * fErrorBandXY -> GetNbinsX()); 00983 00984 if (fFitFunctionIndexX >= 0) 00985 { 00986 for (int ix = 1; ix <= fErrorBandNbinsX; ix++) 00987 { 00988 double sum = 0; 00989 00990 for (int iy = 1; iy <= fErrorBandNbinsY; iy++) 00991 sum += fErrorBandXY -> GetBinContent(ix, iy); 00992 00993 for (int iy = 1; iy <= fErrorBandNbinsY; iy++) 00994 { 00995 double newvalue = fErrorBandXY -> GetBinContent(ix, iy) / sum; 00996 fErrorBandXY -> SetBinContent(ix, iy, newvalue); 00997 } 00998 } 00999 } 01000 01001 BCLog::Out(BCLog::detail, BCLog::detail, 01002 Form("BCIntegrate::MarginalizeAllByMetro done with %i trials and %i accepted trials. Efficiency is %f",fNmetro, fNacceptedMCMC, double(fNacceptedMCMC)/double(fNmetro))); 01003 01004 return fNvar+nh2d; 01005 }
| TH1D * BCIntegrate::GetH1D | ( | int | parIndex | ) |
| parIndex1 | Index of parameter |
Definition at line 1008 of file BCIntegrate.cxx.
01009 { 01010 if(fHProb1D.size()==0) 01011 { 01012 BCLog::Out(BCLog::warning, BCLog::warning, 01013 "BCModel::GetH1D. MarginalizeAll() has to be run prior to this to fill the distributions."); 01014 return 0; 01015 } 01016 01017 if(parIndex<0 || parIndex>fNvar) 01018 { 01019 BCLog::Out(BCLog::warning, BCLog::warning, 01020 Form("BCIntegrate::GetH1D. Parameter index %d is invalid.",parIndex)); 01021 return 0; 01022 } 01023 01024 return fHProb1D[parIndex]; 01025 }
| int BCIntegrate::GetH2DIndex | ( | int | parIndex1, | |
| int | parIndex2 | |||
| ) |
| parIndex1 | Index of first parameter | |
| parIndex2 | Index of second parameter, with parIndex2>parIndex1 |
Definition at line 1028 of file BCIntegrate.cxx.
01029 { 01030 if(parIndex1>fNvar-1 || parIndex1<0) 01031 { 01032 BCLog::Out(BCLog::warning, BCLog::warning, 01033 Form("BCIntegrate::GetH2DIndex. Parameter index %d is invalid", parIndex1)); 01034 return -1; 01035 } 01036 01037 if(parIndex2>fNvar-1 || parIndex2<0) 01038 { 01039 BCLog::Out(BCLog::warning, BCLog::warning, 01040 Form("BCIntegrate::GetH2DIndex. Parameter index %d is invalid", parIndex2)); 01041 return -1; 01042 } 01043 01044 if(parIndex1==parIndex2) 01045 { 01046 BCLog::Out(BCLog::warning, BCLog::warning, 01047 Form("BCIntegrate::GetH2DIndex. Parameters have equal indeces: %d , %d", parIndex1, parIndex2)); 01048 return -1; 01049 } 01050 01051 if(parIndex1>parIndex2) 01052 { 01053 BCLog::Out(BCLog::warning, BCLog::warning, 01054 "BCIntegrate::GetH2DIndex. First parameters must be smaller than second (sorry)."); 01055 return -1; 01056 } 01057 01058 int index=0; 01059 for(int i=0;i<fNvar-1;i++) 01060 for(int j=i+1;j<fNvar;j++) 01061 { 01062 if(i==parIndex1 && j==parIndex2) 01063 return index; 01064 index++; 01065 } 01066 01067 BCLog::Out(BCLog::warning, BCLog::warning, 01068 "BCIntegrate::GetH2DIndex. Invalid index combination."); 01069 01070 return -1; 01071 }
| TH2D * BCIntegrate::GetH2D | ( | int | parIndex1, | |
| int | parIndex2 | |||
| ) |
| parIndex1 | Index of first parameter | |
| parIndex2 | Index of second parameter, with parIndex2>parIndex1 |
Definition at line 1074 of file BCIntegrate.cxx.
01075 { 01076 if(fHProb2D.size()==0) 01077 { 01078 BCLog::Out(BCLog::warning, BCLog::warning, 01079 "BCModel::GetH2D. MarginalizeAll() has to be run prior to this to fill the distributions."); 01080 return 0; 01081 } 01082 01083 int hindex = this -> GetH2DIndex(parIndex1, parIndex2); 01084 if(hindex==-1) 01085 return 0; 01086 01087 if(hindex>(int)fHProb2D.size()-1) 01088 { 01089 BCLog::Out(BCLog::warning, BCLog::warning, 01090 "BCIntegrate::GetH2D. Got invalid index from GetH2DIndex(). Something went wrong."); 01091 return 0; 01092 } 01093 01094 return fHProb2D[hindex]; 01095 }
| void BCIntegrate::InitMetro | ( | ) |
Initializes the Metropolis algorithm (for details see manual)
Definition at line 1130 of file BCIntegrate.cxx.
01131 { 01132 fNmetro=0; 01133 01134 if (fXmetro0.size() <= 0) 01135 { 01136 // start in the center of the phase space 01137 for(int i=0;i<fNvar;i++) 01138 fXmetro0.push_back((fMin[i]+fMax[i])/2.0); 01139 } 01140 01141 fMarkovChainValue = this -> LogEval(fXmetro0); 01142 01143 // run metropolis for a few times and dump the result... (to forget the initial position) 01144 std::vector <double> x; 01145 x.assign(fNvar,0.); 01146 01147 for(int i=0;i<1000;i++) 01148 GetRandomPointMetro(x); 01149 01150 fNmetro = 0; 01151 }
| void BCIntegrate::SAInitialize | ( | ) |
| void BCIntegrate::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 in BCModel.
Definition at line 1304 of file BCIntegrate.cxx.
01305 { 01306 bool have_start = true; 01307 01308 // check start values 01309 if (int(start.size()) != fNvar) 01310 have_start = false; 01311 01312 // set global this 01313 global_this = this; 01314 01315 // define minuit 01316 if (fMinuit) 01317 delete fMinuit; 01318 fMinuit = new TMinuit(fNvar); 01319 01320 // set print level 01321 fMinuit -> SetPrintLevel(printlevel); 01322 01323 // set function 01324 fMinuit -> SetFCN(&BCIntegrate::FCNLikelihood); 01325 01326 // set UP for likelihood 01327 fMinuit -> SetErrorDef(0.5); 01328 01329 // set parameters 01330 int flag; 01331 for (int i = 0; i < fNvar; i++) 01332 { 01333 double starting_point = (fMin[i]+fMax[i])/2.; 01334 if(have_start) 01335 starting_point = start[i]; 01336 fMinuit -> mnparm(i, 01337 fx -> at(i) -> GetName().data(), 01338 starting_point, 01339 (fMax[i]-fMin[i])/100.0, 01340 fMin[i], 01341 fMax[i], 01342 flag); 01343 } 01344 01345 // do mcmc minimization 01346 // fMinuit -> mnseek(); 01347 01348 // do minimization 01349 fMinuit -> mnexcm("MIGRAD", fMinuitArglist, 2, flag); 01350 01351 // improve search for local minimum 01352 // fMinuit -> mnimpr(); 01353 01354 // copy flag 01355 fMinuitErrorFlag = flag; 01356 01357 // check if mode found by minuit is better than previous estimation 01358 double probmax = 0; 01359 bool valid = false; 01360 01361 if ( int(fBestFitParameters.size()) == fNvar) 01362 { 01363 valid = true; 01364 for (int i = 0; i < fNvar; ++i) 01365 if (fBestFitParameters.at(i) < fMin[i] || fBestFitParameters.at(i) > fMax[i]) 01366 valid= false; 01367 01368 if (valid) 01369 probmax = this -> Eval( fBestFitParameters ); 01370 } 01371 01372 std::vector<double> tempvec; 01373 for (int i = 0; i < fNvar; i++) 01374 { 01375 double par; 01376 double parerr; 01377 fMinuit -> GetParameter(i, par, parerr); 01378 tempvec.push_back(par); 01379 } 01380 double probmaxminuit = this -> Eval( tempvec ); 01381 01382 // set best fit parameters 01383 if ( (probmaxminuit > probmax && !fFlagIgnorePrevOptimization) || !valid || fFlagIgnorePrevOptimization) 01384 { 01385 fBestFitParameters.clear(); 01386 fBestFitParameterErrors.clear(); 01387 01388 for (int i = 0; i < fNvar; i++) 01389 { 01390 double par; 01391 double parerr; 01392 fMinuit -> GetParameter(i, par, parerr); 01393 fBestFitParameters.push_back(par); 01394 fBestFitParameterErrors.push_back(parerr); 01395 } 01396 01397 // set optimization method used to find the mode 01398 fOptimizationMethodMode = BCIntegrate::kOptMinuit; 01399 } 01400 01401 // delete minuit 01402 delete fMinuit; 01403 fMinuit = 0; 01404 01405 return; 01406 }
| void BCIntegrate::FindModeMCMC | ( | ) |
Does the mode finding using Markov Chain Monte Carlo
Definition at line 1821 of file BCIntegrate.cxx.
01822 { 01823 // call PreRun 01824 // if (flag_run == 0) 01825 // if (!fMCMCFlagPreRun) 01826 // this -> MCMCMetropolisPreRun(); 01827 01828 // find global maximum 01829 // double probmax = (this -> MCMCGetMaximumLogProb()).at(0); 01830 01831 double probmax = 0; 01832 01833 if ( int(fBestFitParameters.size()) == fNvar) 01834 { 01835 probmax = this -> Eval( fBestFitParameters ); 01836 // fBestFitParameters = this -> MCMCGetMaximumPoint(0); 01837 } 01838 01839 // loop over all chains and find the maximum point 01840 for (int i = 0; i < fMCMCNChains; ++i) 01841 { 01842 double prob = exp( (this -> MCMCGetMaximumLogProb()).at(i)); 01843 01844 // copy the point into the vector 01845 if ( (prob >= probmax && !fFlagIgnorePrevOptimization) || fFlagIgnorePrevOptimization) 01846 { 01847 probmax = prob; 01848 01849 fBestFitParameters.clear(); 01850 fBestFitParameterErrors.clear(); 01851 fBestFitParameters = this -> MCMCGetMaximumPoint(i); 01852 01853 for (int j = 0; j < fNvar; ++j) 01854 fBestFitParameterErrors.push_back(0.0); 01855 01856 fOptimizationMethodMode = BCIntegrate::kOptMetropolis; 01857 } 01858 } 01859 01860 }
| void BCIntegrate::FindModeSA | ( | std::vector< double > | start = std::vector< double >(0) |
) |
Does the mode finding using Simulated Annealing. If starting point is not specified, finding will start from the center of the parameter space.
| start | point in parameter space from thich the mode finding is started. |
Definition at line 1410 of file BCIntegrate.cxx.
01411 { 01412 // note: if f(x) is the function to be minimized, then 01413 // f(x) := - this->LogEval(parameters) 01414 01415 bool have_start = true; 01416 std::vector<double> x, y, best_fit; // vectors for current state, new proposed state and best fit up to now 01417 double fval_x, fval_y, fval_best_fit; // function values at points x, y and best_fit (we save them rather than to re-calculate them every time) 01418 int t = 1; // time iterator 01419 01420 // check start values 01421 if (int(start.size()) != fNvar) 01422 have_start = false; 01423 01424 // if no starting point is given, set to center of parameter space 01425 if ( !have_start ) 01426 { 01427 start.clear(); 01428 for (int i = 0; i < fNvar; i++) 01429 start.push_back((fMin[i]+fMax[i])/2.); 01430 } 01431 01432 // set current state and best fit to starting point 01433 x.clear(); 01434 best_fit.clear(); 01435 for (int i = 0; i < fNvar; i++) 01436 { 01437 x.push_back(start[i]); 01438 best_fit.push_back(start[i]); 01439 } 01440 // calculate function value at starting point 01441 fval_x = fval_best_fit = this->LogEval(x); 01442 01443 // run while still "hot enough" 01444 while ( this->SATemperature(t) > fSATmin ) 01445 { 01446 // generate new state 01447 y = this->GetProposalPointSA(x, t); 01448 01449 // check if the proposed point is inside the phase space 01450 // if not, reject it 01451 bool is_in_ranges = true; 01452 01453 for (int i = 0; i < fNvar; i++) 01454 if (y[i] > fMax[i] || y[i] < fMin[i]) 01455 is_in_ranges = false; 01456 01457 if ( !is_in_ranges ) 01458 ; // do nothing... 01459 else 01460 { 01461 // calculate function value at new point 01462 fval_y = this->LogEval(y); 01463 01464 // is it better than the last one? 01465 // if so, update state and chef if it is the new best fit... 01466 if (fval_y >= fval_x) 01467 { 01468 x.clear(); 01469 for (int i = 0; i < fNvar; i++) 01470 x.push_back(y[i]); 01471 01472 fval_x = fval_y; 01473 01474 if (fval_y > fval_best_fit) 01475 { 01476 best_fit.clear(); 01477 for (int i = 0; i < fNvar; i++) 01478 best_fit.push_back(y[i]); 01479 01480 fval_best_fit = fval_y; 01481 } 01482 } 01483 // ...else, only accept new state w/ certain probability 01484 else 01485 { 01486 if (fRandom->Rndm() <= exp( (fval_y - fval_x) 01487 / this->SATemperature(t) )) 01488 { 01489 x.clear(); 01490 for (int i = 0; i < fNvar; i++) 01491 x.push_back(y[i]); 01492 01493 fval_x = fval_y; 01494 } 01495 } 01496 } 01497 01498 // update tree variables 01499 fSANIterations = t; 01500 fSATemperature = this -> SATemperature(t); 01501 fSALogProb = fval_x; 01502 fSAx.clear(); 01503 for (int i = 0; i < fNvar; ++i) 01504 fSAx.push_back(x[i]); 01505 01506 // fill tree 01507 if (fFlagWriteSAToFile) 01508 fTreeSA -> Fill(); 01509 01510 // increate t 01511 t++; 01512 } 01513 01514 // check if mode found by minuit is better than previous estimation 01515 double probmax = 0; 01516 bool valid = false; 01517 01518 if ( int(fBestFitParameters.size()) == fNvar) 01519 { 01520 valid = true; 01521 for (int i = 0; i < fNvar; ++i) 01522 if (fBestFitParameters.at(i) < fMin[i] || fBestFitParameters.at(i) > fMax[i]) 01523 valid= false; 01524 01525 if (valid) 01526 probmax = this -> Eval( fBestFitParameters ); 01527 } 01528 01529 double probmaxsa = this -> Eval( best_fit); 01530 01531 if ( (probmaxsa > probmax && !fFlagIgnorePrevOptimization) || !valid || fFlagIgnorePrevOptimization) 01532 { 01533 // set best fit parameters 01534 fBestFitParameters.clear(); 01535 fBestFitParameterErrors.clear(); 01536 01537 for (int i = 0; i < fNvar; i++) 01538 { 01539 fBestFitParameters.push_back(best_fit[i]); 01540 fBestFitParameterErrors.push_back(0.0); 01541 } 01542 01543 // set optimization moethod used to find the mode 01544 fOptimizationMethodMode = BCIntegrate::kOptSA; 01545 } 01546 01547 return; 01548 }
| double BCIntegrate::SATemperature | ( | double | t | ) |
Temperature annealing schedule for use with Simulated Annealing. Delegates to the appropriate method according to fSASchedule.
| t | iterator for lowering the temperature over time. |
Definition at line 1552 of file BCIntegrate.cxx.
01553 { 01554 // do we have Cauchy (default), Boltzmann or custom annealing schedule? 01555 if (this->fSASchedule == BCIntegrate::kSABoltzmann) 01556 return this->SATemperatureBoltzmann(t); 01557 else if (this->fSASchedule == BCIntegrate::kSACauchy) 01558 return this->SATemperatureCauchy(t); 01559 else 01560 return this->SATemperatureCustom(t); 01561 }
| double BCIntegrate::SATemperatureBoltzmann | ( | double | t | ) |
Temperature annealing schedule for use with Simulated Annealing. This method is used for Boltzmann annealing schedule.
| t | iterator for lowering the temperature over time. |
Definition at line 1565 of file BCIntegrate.cxx.
01566 { 01567 return fSAT0 / log((double)(t + 1)); 01568 }
| double BCIntegrate::SATemperatureCauchy | ( | double | t | ) |
Temperature annealing schedule for use with Simulated Annealing. This method is used for Cauchy annealing schedule.
| t | iterator for lowering the temperature over time. |
Definition at line 1572 of file BCIntegrate.cxx.
01573 { 01574 return fSAT0 / (double)t; 01575 }
| double BCIntegrate::SATemperatureCustom | ( | double | t | ) | [virtual] |
Temperature annealing schedule for use with Simulated Annealing. This is a virtual method to be overridden by a user-defined custom temperature schedule.
| t | iterator for lowering the temperature over time. |
Definition at line 1579 of file BCIntegrate.cxx.
01580 { 01581 BCLog::OutError("BCIntegrate::SATemperatureCustom : No custom temperature schedule defined"); 01582 return 0.; 01583 }
| std::vector< double > BCIntegrate::GetProposalPointSA | ( | std::vector< double > | x, | |
| int | t | |||
| ) |
Generates a new state in a neighbourhood around x that is to be accepted or rejected by the Simulated Annealing algorithm. Delegates the generation to the appropriate method according to fSASchedule.
| x | last state. | |
| t | time iterator to determine current temperature. |
Definition at line 1587 of file BCIntegrate.cxx.
01588 { 01589 // do we have Cauchy (default), Boltzmann or custom annealing schedule? 01590 if (this->fSASchedule == BCIntegrate::kSABoltzmann) 01591 return this->GetProposalPointSABoltzmann(x, t); 01592 else if (this->fSASchedule == BCIntegrate::kSACauchy) 01593 return this->GetProposalPointSACauchy(x, t); 01594 else 01595 return this->GetProposalPointSACustom(x, t); 01596 }
| std::vector< double > BCIntegrate::GetProposalPointSABoltzmann | ( | std::vector< double > | x, | |
| int | t | |||
| ) |
Generates a new state in a neighbourhood around x that is to be accepted or rejected by the Simulated Annealing algorithm. This method is used for Boltzmann annealing schedule.
| x | last state. | |
| t | time iterator to determine current temperature. |
Definition at line 1600 of file BCIntegrate.cxx.
01601 { 01602 std::vector<double> y; 01603 y.clear(); 01604 double new_val, norm; 01605 01606 for (int i = 0; i < fNvar; i++) 01607 { 01608 norm = (fMax[i] - fMin[i]) * this->SATemperature(t) / 2.; 01609 new_val = x[i] + norm * fRandom->Gaus(); 01610 y.push_back(new_val); 01611 } 01612 01613 return y; 01614 }
| std::vector< double > BCIntegrate::GetProposalPointSACauchy | ( | std::vector< double > | x, | |
| int | t | |||
| ) |
Generates a new state in a neighbourhood around x that is to be accepted or rejected by the Simulated Annealing algorithm. This method is used for Cauchy annealing schedule.
| x | last state. | |
| t | time iterator to determine current temperature. |
Definition at line 1618 of file BCIntegrate.cxx.
01619 { 01620 std::vector<double> y; 01621 y.clear(); 01622 01623 if (fNvar == 1) 01624 { 01625 double cauchy, new_val, norm; 01626 01627 norm = (fMax[0] - fMin[0]) * this->SATemperature(t) / 2.; 01628 cauchy = tan(3.14159 * (fRandom->Rndm() - 0.5)); 01629 new_val = x[0] + norm * cauchy; 01630 y.push_back(new_val); 01631 } 01632 else 01633 { 01634 // use sampling to get radial n-dim Cauchy distribution 01635 01636 // first generate a random point uniformly distributed on a 01637 // fNvar-dimensional hypersphere 01638 y = this->SAHelperGetRandomPointOnHypersphere(); 01639 01640 // scale the vector by a random factor determined by the radial 01641 // part of the fNvar-dimensional Cauchy distribution 01642 double radial = this->SATemperature(t) * this->SAHelperGetRadialCauchy(); 01643 01644 // scale y by radial part and the size of dimension i in phase space 01645 // afterwards, move by x 01646 for (int i = 0; i < fNvar; i++) 01647 y[i] = (fMax[i] - fMin[i]) * y[i] * radial / 2. + x[i]; 01648 } 01649 01650 return y; 01651 }
| std::vector< double > BCIntegrate::GetProposalPointSACustom | ( | std::vector< double > | x, | |
| int | t | |||
| ) | [virtual] |
Generates a new state in a neighbourhood around x that is to be accepted or rejected by the Simulated Annealing algorithm. This is a virtual method to be overridden by a user-defined custom proposal function.
| x | last state. | |
| t | time iterator to determine current temperature. |
Definition at line 1655 of file BCIntegrate.cxx.
01656 { 01657 BCLog::OutError("BCIntegrate::GetProposalPointSACustom : No custom proposal function defined"); 01658 return std::vector<double>(fNvar); 01659 }
| std::vector< double > BCIntegrate::SAHelperGetRandomPointOnHypersphere | ( | ) |
Generates a uniform distributed random point on the surface of a fNvar-dimensional Hypersphere. Used as a helper to generate proposal points for Cauchy annealing.
Definition at line 1663 of file BCIntegrate.cxx.
01664 { 01665 std::vector<double> rand_point(fNvar); 01666 01667 // This method can only be called with fNvar >= 2 since the 1-dim case 01668 // is already hard wired into the Cauchy annealing proposal function. 01669 // To speed things up, hard-code fast method for 2 and dimensions. 01670 // The algorithm for 2D can be found at 01671 // http://mathworld.wolfram.com/CirclePointPicking.html 01672 // For 3D just using ROOT's algorithm. 01673 if (fNvar == 2) 01674 { 01675 double x1, x2, s; 01676 do { 01677 x1 = fRandom->Rndm() * 2. - 1.; 01678 x2 = fRandom->Rndm() * 2. - 1.; 01679 s = x1*x1 + x2*x2; 01680 } while (s >= 1); 01681 01682 rand_point[0] = (x1*x1 - x2*x2) / s; 01683 rand_point[1] = (2.*x1*x2) / s; 01684 } 01685 else if (fNvar == 3) 01686 { 01687 fRandom->Sphere(rand_point[0], rand_point[1], rand_point[2], 1.0); 01688 } 01689 else 01690 { 01691 double s = 0., 01692 gauss_num; 01693 01694 for (int i = 0; i < fNvar; i++) 01695 { 01696 gauss_num = fRandom->Gaus(); 01697 rand_point[i] = gauss_num; 01698 s += gauss_num * gauss_num; 01699 } 01700 s = sqrt(s); 01701 01702 for (int i = 0; i < fNvar; i++) 01703 rand_point[i] = rand_point[i] / s; 01704 } 01705 01706 return rand_point; 01707 }
| double BCIntegrate::SAHelperGetRadialCauchy | ( | ) |
Generates the radial part of a n-dimensional Cauchy distribution. Helper function for Cauchy annealing.
Definition at line 1711 of file BCIntegrate.cxx.
01712 { 01713 // theta is sampled from a rather complicated distribution, 01714 // so first we create a lookup table with 10000 random numbers 01715 // once and then, each time we need a new random number, 01716 // we just look it up in the table. 01717 double theta; 01718 01719 // static vectors for theta-sampling-map 01720 static double map_u[10001]; 01721 static double map_theta[10001]; 01722 static bool initialized = false; 01723 static int map_dimension = 0; 01724 01725 // is the lookup-table already initialized? if not, do it! 01726 if (!initialized || map_dimension != fNvar) 01727 { 01728 double init_theta; 01729 double beta = this->SAHelperSinusToNIntegral(fNvar - 1, 1.57079632679); 01730 01731 for (int i = 0; i <= 10000; i++) 01732 { 01733 init_theta = 3.14159265 * (double)i / 5000.; 01734 map_theta[i] = init_theta; 01735 01736 map_u[i] = this->SAHelperSinusToNIntegral(fNvar - 1, init_theta) / beta; 01737 } 01738 01739 map_dimension = fNvar; 01740 initialized = true; 01741 } // initializing is done. 01742 01743 // generate uniform random number for sampling 01744 double u = fRandom->Rndm(); 01745 01746 // Find the two elements just greater than and less than u 01747 // using a binary search (O(log(N))). 01748 int lo = 0; 01749 int up = 10000; 01750 int mid; 01751 01752 while (up != lo) 01753 { 01754 mid = ((up - lo + 1) / 2) + lo; 01755 01756 if (u >= map_u[mid]) 01757 lo = mid; 01758 else 01759 up = mid - 1; 01760 } 01761 up++; 01762 01763 // perform linear interpolation: 01764 theta = map_theta[lo] + (u - map_u[lo]) / (map_u[up] - map_u[lo]) 01765 * (map_theta[up] - map_theta[lo]); 01766 01767 return tan(theta); 01768 }
| double BCIntegrate::SAHelperSinusToNIntegral | ( | int | dim, | |
| double | theta | |||
| ) |
Returns the Integral of sin^dim from 0 to theta. Helper function needed for generating Cauchy distributions.
Definition at line 1772 of file BCIntegrate.cxx.
01773 { 01774 if (dim < 1) 01775 return theta; 01776 else if (dim == 1) 01777 return (1. - cos(theta)); 01778 else if (dim == 2) 01779 return 0.5 * (theta - sin(theta) * cos(theta)); 01780 else if (dim == 3) 01781 return (2. - sin(theta) * sin(theta) * cos(theta) - 2. * cos(theta)) / 3.; 01782 else 01783 return - pow(sin(theta), (double)(dim - 1)) * cos(theta) / (double)dim 01784 + (double)(dim - 1) / (double)dim 01785 * this->SAHelperSinusToNIntegral(dim - 2, theta); 01786 }
| void BCIntegrate::FCNLikelihood | ( | int & | npar, | |
| double * | grad, | |||
| double & | fval, | |||
| double * | par, | |||
| int | flag | |||
| ) | [static] |
Definition at line 1802 of file BCIntegrate.cxx.
01803 { 01804 // copy parameters 01805 std::vector <double> parameters; 01806 01807 int n = global_this -> GetNvar(); 01808 01809 for (int i = 0; i < n; i++) 01810 parameters.push_back(par[i]); 01811 01812 fval = - global_this -> LogEval(parameters); 01813 01814 // delete parameters 01815 parameters.clear(); 01816 }
| virtual void BCIntegrate::MCMCUserIterationInterface | ( | ) | [inline, virtual] |
| void BCIntegrate::MCMCIterationInterface | ( | ) | [private, virtual] |
Reimplemented from BCEngineMCMC.
Definition at line 2043 of file BCIntegrate.cxx.
02044 { 02045 // what's within this method will be executed 02046 // for every iteration of the MCMC 02047 02048 // fill error band 02049 this -> MCMCFillErrorBand(); 02050 02051 // do user defined stuff 02052 this -> MCMCUserIterationInterface(); 02053 }
| void BCIntegrate::MCMCFillErrorBand | ( | ) | [private] |
Definition at line 2056 of file BCIntegrate.cxx.
02057 { 02058 if (!fFillErrorBand) 02059 return; 02060 02061 // function fitting 02062 if (fFitFunctionIndexX < 0) 02063 return; 02064 02065 // loop over all possible x values ... 02066 if (fErrorBandContinuous) 02067 { 02068 double x = 0; 02069 for (int ix = 0; ix < fErrorBandNbinsX; ix++) 02070 { 02071 // calculate x 02072 x = fErrorBandXY -> GetXaxis() -> GetBinCenter(ix + 1); 02073 02074 // calculate y 02075 std::vector <double> xvec; 02076 xvec.push_back(x); 02077 02078 // loop over all chains 02079 for (int ichain = 0; ichain < this -> MCMCGetNChains(); ++ichain) 02080 { 02081 // calculate y 02082 double y = this -> FitFunction(xvec, this -> MCMCGetx(ichain)); 02083 02084 // fill histogram 02085 fErrorBandXY -> Fill(x, y); 02086 } 02087 02088 xvec.clear(); 02089 } 02090 } 02091 // ... or evaluate at the data point x-values 02092 else 02093 { 02094 int ndatapoints = int(fErrorBandX.size()); 02095 double x = 0; 02096 02097 for (int ix = 0; ix < ndatapoints; ++ix) 02098 { 02099 // calculate x 02100 x = fErrorBandX.at(ix); 02101 02102 // calculate y 02103 std::vector <double> xvec; 02104 xvec.push_back(x); 02105 02106 // loop over all chains 02107 for (int ichain = 0; ichain < this -> MCMCGetNChains(); ++ichain) 02108 { 02109 // calculate y 02110 double y = this -> FitFunction(xvec, this -> MCMCGetx(ichain)); 02111 02112 // fill histogram 02113 fErrorBandXY -> Fill(x, y); 02114 } 02115 02116 xvec.clear(); 02117 } 02118 } 02119 }
BCParameterSet* BCIntegrate::fx [private] |
Set of parameters for the integration.
Definition at line 719 of file BCIntegrate.h.
double* BCIntegrate::fMin [private] |
Array containing the lower boundaries of the variables to integrate over.
Definition at line 731 of file BCIntegrate.h.
double* BCIntegrate::fMax [private] |
Array containing the upper boundaries of the variables to integrate over.
Definition at line 735 of file BCIntegrate.h.
int* BCIntegrate::fVarlist [private] |
List of variables containing a flag whether to integrate over them or not.
Definition at line 739 of file BCIntegrate.h.
int BCIntegrate::fNiterPerDimension [private] |
Number of iteration per dimension for Monte Carlo integration.
Definition at line 743 of file BCIntegrate.h.
Current integration method
Definition at line 747 of file BCIntegrate.h.
Current marginalization method
Definition at line 751 of file BCIntegrate.h.
Current mode finding method
Definition at line 755 of file BCIntegrate.h.
Method with which the global mode was found (can differ from fOptimization method in case more than one algorithm is used).
Definition at line 760 of file BCIntegrate.h.
Current Simulated Annealing schedule
Definition at line 764 of file BCIntegrate.h.
int BCIntegrate::fNIterationsMax [private] |
Maximum number of iterations
Definition at line 768 of file BCIntegrate.h.
int BCIntegrate::fNIterations [private] |
Number of iterations in the most recent Monte Carlo integation
Definition at line 772 of file BCIntegrate.h.
double BCIntegrate::fRelativePrecision [private] |
Relative precision aimed at in the Monte Carlo integation
Definition at line 776 of file BCIntegrate.h.
double BCIntegrate::fError [private] |
The uncertainty in the most recent Monte Carlo integration
Definition at line 780 of file BCIntegrate.h.
int BCIntegrate::fNmetro [private] |
The number of iterations in the Metropolis integration
Definition at line 784 of file BCIntegrate.h.
int BCIntegrate::fNacceptedMCMC [private] |
Definition at line 785 of file BCIntegrate.h.
std::vector<double> BCIntegrate::fXmetro0 [private] |
A vector of points in parameter space used for the Metropolis algorithm
Definition at line 789 of file BCIntegrate.h.
std::vector<double> BCIntegrate::fXmetro1 [private] |
A vector of points in parameter space used for the Metropolis algorithm
Definition at line 793 of file BCIntegrate.h.
double BCIntegrate::fMarkovChainValue [private] |
A double containing the log likelihood value at the point fXmetro1
Definition at line 797 of file BCIntegrate.h.
int BCIntegrate::fNvar [protected] |
Number of variables to integrate over.
Definition at line 811 of file BCIntegrate.h.
int BCIntegrate::fNbins [protected] |
Number of bins per dimension for the marginalized distributions
Definition at line 815 of file BCIntegrate.h.
int BCIntegrate::fNSamplesPer2DBin [protected] |
Number of samples per 2D bin per variable in the Metropolis marginalization.
Definition at line 820 of file BCIntegrate.h.
double BCIntegrate::fMarkovChainStepSize [protected] |
Step size in the Markov chain relative to min and max
Definition at line 824 of file BCIntegrate.h.
int BCIntegrate::fMarkovChainNIterations [protected] |
Definition at line 826 of file BCIntegrate.h.
bool BCIntegrate::fMarkovChainAutoN [protected] |
Definition at line 828 of file BCIntegrate.h.
BCDataPoint* BCIntegrate::fDataPointLowerBoundaries [protected] |
data point containing the lower boundaries of possible data values
Definition at line 832 of file BCIntegrate.h.
BCDataPoint* BCIntegrate::fDataPointUpperBoundaries [protected] |
data point containing the upper boundaries of possible data values
Definition at line 836 of file BCIntegrate.h.
std::vector<bool> BCIntegrate::fDataFixedValues [protected] |
Definition at line 838 of file BCIntegrate.h.
TRandom3* BCIntegrate::fRandom [protected] |
A ROOT random number generator
Definition at line 842 of file BCIntegrate.h.
std::vector<double> BCIntegrate::fBestFitParameters [protected] |
A vector of best fit parameters estimated from the global probability and the estimate on their uncertainties
Definition at line 847 of file BCIntegrate.h.
std::vector<double> BCIntegrate::fBestFitParameterErrors [protected] |
Definition at line 848 of file BCIntegrate.h.
std::vector<double> BCIntegrate::fBestFitParametersMarginalized [protected] |
A vector of best fit parameters estimated from the marginalized probability
Definition at line 852 of file BCIntegrate.h.
std::vector<TH1D *> BCIntegrate::fHProb1D [protected] |
Vector of TH1D histograms for marginalized probability distributions
Definition at line 856 of file BCIntegrate.h.
std::vector<TH2D *> BCIntegrate::fHProb2D [protected] |
Vector of TH2D histograms for marginalized probability distributions
Definition at line 860 of file BCIntegrate.h.
bool BCIntegrate::fFillErrorBand [protected] |
Flag whether or not to fill the error band
Definition at line 864 of file BCIntegrate.h.
int BCIntegrate::fFitFunctionIndexX [protected] |
The indeces for function fits
Definition at line 868 of file BCIntegrate.h.
int BCIntegrate::fFitFunctionIndexY [protected] |
Definition at line 869 of file BCIntegrate.h.
bool BCIntegrate::fErrorBandContinuous [protected] |
A flag for single point evaluation of the error "band"
Definition at line 873 of file BCIntegrate.h.
std::vector<double> BCIntegrate::fErrorBandX [protected] |
Definition at line 874 of file BCIntegrate.h.
TH2D* BCIntegrate::fErrorBandXY [protected] |
The error band histogram
Definition at line 878 of file BCIntegrate.h.
int BCIntegrate::fErrorBandNbinsX [protected] |
Number of X bins of the error band histogram
Definition at line 882 of file BCIntegrate.h.
int BCIntegrate::fErrorBandNbinsY [protected] |
Nnumber of Y bins of the error band histogram
Definition at line 886 of file BCIntegrate.h.
TMinuit* BCIntegrate::fMinuit [protected] |
Minuit
Definition at line 890 of file BCIntegrate.h.
double BCIntegrate::fMinuitArglist[2] [protected] |
Definition at line 892 of file BCIntegrate.h.
int BCIntegrate::fMinuitErrorFlag [protected] |
Definition at line 893 of file BCIntegrate.h.
double BCIntegrate::fFlagIgnorePrevOptimization [protected] |
Flag for ignoring older results of minimization
Definition at line 897 of file BCIntegrate.h.
bool BCIntegrate::fFlagWriteMarkovChain [protected] |
Flag for writing Markov chain to file
Definition at line 901 of file BCIntegrate.h.
TTree* BCIntegrate::fMarkovChainTree [protected] |
ROOT tree containing the Markov chain
Definition at line 905 of file BCIntegrate.h.
int BCIntegrate::fMCMCIteration [protected] |
Iteration of the MCMC
Definition at line 909 of file BCIntegrate.h.
double BCIntegrate::fSAT0 [protected] |
Starting temperature for Simulated Annealing
Definition at line 913 of file BCIntegrate.h.
double BCIntegrate::fSATmin [protected] |
Minimal/Threshold temperature for Simulated Annealing
Definition at line 917 of file BCIntegrate.h.
TTree* BCIntegrate::fTreeSA [protected] |
Tree for the Simulated Annealing
Definition at line 921 of file BCIntegrate.h.
bool BCIntegrate::fFlagWriteSAToFile [protected] |
Flag deciding whether to write SA to file or not.
Definition at line 925 of file BCIntegrate.h.
int BCIntegrate::fSANIterations [protected] |
Definition at line 927 of file BCIntegrate.h.
double BCIntegrate::fSATemperature [protected] |
Definition at line 928 of file BCIntegrate.h.
double BCIntegrate::fSALogProb [protected] |
Definition at line 929 of file BCIntegrate.h.
std::vector<double> BCIntegrate::fSAx [protected] |
Definition at line 930 of file BCIntegrate.h.
1.5.1