BCHistogramFitter Class Reference

A class for fitting histograms with functions. More...

#include <BCHistogramFitter.h>

Inherits BCModel.

Collaboration diagram for BCHistogramFitter:
Collaboration graph
[legend]

List of all members.

Public Member Functions

Constructors and destructors



 BCHistogramFitter (TH1D *hist, TF1 *func)
 BCHistogramFitter ()
 ~BCHistogramFitter ()
Member functions (miscellaneous methods)



int CalculatePValueFast (std::vector< double > par, double &pvalue)
void DrawFit (const char *options="", bool flaglegend=false)
int Fit (TH1D *hist, TF1 *func)
int Fit ()
double FitFunction (std::vector< double > x, std::vector< double > parameters)
virtual double LogAPrioriProbability (std::vector< double > parameters)
virtual double LogLikelihood (std::vector< double > parameters)
Member functions (get)



TGraph * GetErrorBand ()
TF1 * GetFitFunction ()
TGraph * GetGraphFitFunction ()
TH1D * GetHistogram ()
Member functions (set)



int SetFitFunction (TF1 *func)
void SetFlagIntegration (bool flag)
int SetHistogram (TH1D *hist)

Private Attributes

TGraph * fErrorBand
TF1 * fFitFunction
bool fFlagIntegration
TGraph * fGraphFitFunction
TH1D * fHistogram

Detailed Description

A class for fitting histograms with functions.

Author:
Daniel Kollar
Kevin Kröninger
Version:
1.0
Date:
11.2008 This class allows fitting of a TH1D histogram using a TF1 function.

Definition at line 34 of file BCHistogramFitter.h.

Constructor & Destructor Documentation

BCHistogramFitter::BCHistogramFitter (  ) 

The default constructor.

Definition at line 27 of file BCHistogramFitter.cxx.

00027                                      : BCModel("HistogramFitter")
00028 {
00029    fHistogram = 0;
00030    fFitFunction = 0;
00031 
00032    // set default options and values
00033    this -> MCMCSetNIterationsRun(2000);
00034    this -> SetFillErrorBand();
00035    fFlagIntegration = true;
00036 }

BCHistogramFitter::BCHistogramFitter ( TH1D *  hist,
TF1 *  func 
)

A constructor.

Parameters:
hist The histogram (TH1D).
func The fit function.

Definition at line 40 of file BCHistogramFitter.cxx.

00040                                                             : BCModel("HistogramFitter")
00041 {
00042    fHistogram = 0;
00043    fFitFunction = 0;
00044 
00045    this -> SetHistogram(hist);
00046    this -> SetFitFunction(func);
00047 
00048    this -> MCMCSetNIterationsRun(2000);
00049    this -> SetFillErrorBand();
00050    fFlagIntegration = true;
00051 }

BCHistogramFitter::~BCHistogramFitter (  ) 

The default destructor.

Definition at line 147 of file BCHistogramFitter.cxx.

00148 {}

Member Function Documentation

int BCHistogramFitter::CalculatePValueFast ( std::vector< double >  par,
double &  pvalue 
)

Calculate the p-value using fast-MCMC.

Parameters:
par A set of parameter values
pvalue The pvalue
Returns:
An error code

Definition at line 334 of file BCHistogramFitter.cxx.

00335 {
00336    // check size of parameter vector
00337    if (par.size() != this -> GetNParameters())
00338    {
00339       BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::CalculatePValueFast() : Number of parameters is inconsistent.");
00340       return 0;
00341    }
00342 
00343    // check if histogram exists
00344    if (!fHistogram)
00345    {
00346       BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::CalculatePValueFast() : Histogram not defined.");
00347       return 0;
00348    }
00349 
00350    // define temporary variables
00351    int nbins = fHistogram -> GetNbinsX();
00352 
00353    std::vector <int> histogram;
00354    std::vector <double> expectation;
00355    histogram.assign(nbins, 0);
00356    expectation.assign(nbins, 0);
00357 
00358    double logp = 0;
00359    double logp_start = 0;
00360    int counter_pvalue = 0;
00361 
00362    // define starting distribution
00363    for (int ibin = 0; ibin < nbins; ++ibin)
00364    {
00365       // get bin boundaries
00366       double xmin = fHistogram -> GetBinLowEdge(ibin+1);
00367       double xmax = fHistogram -> GetBinLowEdge(ibin+2);
00368 
00369       // get the number of expected events
00370       double yexp = fFitFunction -> Integral(xmin, xmax);
00371 
00372       histogram[ibin]   = int(yexp);
00373       expectation[ibin] = yexp;
00374 
00375       // calculate p;
00376       logp += BCMath::LogPoisson(double(int(yexp)), yexp);
00377       logp_start += BCMath::LogPoisson(fHistogram -> GetBinContent(ibin+1), yexp);
00378    }
00379 
00380    int niter = 100000;
00381 
00382    // loop over iterations
00383    for (int iiter = 0; iiter < niter; ++iiter)
00384    {
00385       // loop over bins
00386       for (int ibin = 0; ibin < nbins; ++ibin)
00387       {
00388          // random step up or down in statistics for this bin
00389          double ptest = fRandom -> Rndm() - 0.5;
00390 
00391          // increase statistics by 1
00392          if (ptest > 0)
00393          {
00394             // calculate factor of probability
00395             double r = expectation.at(ibin) / double(histogram.at(ibin) + 1);
00396 
00397             // walk, or don't (this is the Metropolis part)
00398             if (fRandom -> Rndm() < r)
00399             {
00400                histogram[ibin] = histogram.at(ibin) + 1;
00401                logp += log(r);
00402             }
00403          }
00404 
00405          // decrease statistics by 1
00406          else if (ptest <= 0 && histogram[ibin] > 0)
00407          {
00408             // calculate factor of probability
00409             double r = double(histogram.at(ibin)) / expectation.at(ibin);
00410 
00411             // walk, or don't (this is the Metropolis part)
00412             if (fRandom -> Rndm() < r)
00413             {
00414                histogram[ibin] = histogram.at(ibin) - 1;
00415                logp += log(r);
00416             }
00417          }
00418          // debugKK
00419 //       std::cout << " negative " << std::endl;
00420 
00421       } // end of looping over bins
00422 
00423       // increase counter
00424       if (logp < logp_start)
00425          counter_pvalue++;
00426 
00427    } // end of looping over iterations
00428 
00429    // calculate p-value
00430    pvalue = double(counter_pvalue) / double(niter);
00431 
00432    // no error
00433    return 1;
00434 }

void BCHistogramFitter::DrawFit ( const char *  options = "",
bool  flaglegend = false 
)

Draw the fit in the current pad.

Definition at line 283 of file BCHistogramFitter.cxx.

00284 {
00285    if (!fHistogram)
00286    {
00287       BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::DrawFit() : Histogram not defined.");
00288       return;
00289    }
00290 
00291    if (!fFitFunction)
00292    {
00293       BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::DrawFit() : Fit function not defined.");
00294       return;
00295    }
00296 
00297    // check wheather options contain "same"
00298    TString opt = options;
00299    opt.ToLower();
00300 
00301    // if not same, draw the histogram first to get the axes
00302    if(!opt.Contains("same"))
00303       fHistogram -> Draw(opt.Data());
00304 
00305    // draw the error band as central 68% probability interval
00306    fErrorBand = this -> GetErrorBandGraph(0.16, 0.84);
00307    fErrorBand -> Draw("f same");
00308 
00309    // now draw the histogram again since it was covered by the band
00310    fHistogram -> Draw(TString::Format("%ssame",opt.Data()));
00311 
00312    // draw the fit function on top
00313    fGraphFitFunction = this -> GetFitFunctionGraph( this->GetBestFitParameters() );
00314    fGraphFitFunction -> SetLineColor(kRed);
00315    fGraphFitFunction -> SetLineWidth(2);
00316    fGraphFitFunction -> Draw("l same");
00317 
00318    // draw legend
00319    if (flaglegend)
00320    {
00321       TLegend * legend = new TLegend(0.25, 0.75, 0.55, 0.95);
00322       legend -> SetBorderSize(0);
00323       legend -> SetFillColor(kWhite);
00324       legend -> AddEntry(fHistogram, "Data", "L");
00325       legend -> AddEntry(fGraphFitFunction, "Best fit", "L");
00326       legend -> AddEntry(fErrorBand, "Error band", "F");
00327       legend -> Draw();
00328    }
00329 
00330    gPad -> RedrawAxis();
00331 }

int BCHistogramFitter::Fit ( TH1D *  hist,
TF1 *  func 
)

Performs the fit.

Parameters:
hist The histogram (TH1D).
func The fit function.
Returns:
An error code.

Definition at line 240 of file BCHistogramFitter.cxx.

00241 {
00242    // set histogram
00243    if (hist)
00244       this -> SetHistogram(hist);
00245    else
00246    {
00247       BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::Fit() : Histogram not defined.");
00248       return 0;
00249    }
00250 
00251    // set function
00252    if (func)
00253       this -> SetFitFunction(func);
00254    else
00255    {
00256       BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::Fit() : Fit function not defined.");
00257       return 0;
00258    }
00259 
00260    // perform marginalization
00261    this -> MarginalizeAll();
00262 
00263    // maximize posterior probability, using the best-fit values close
00264    // to the global maximum from the MCMC
00265    this -> FindModeMinuit( this -> GetBestFitParameters() , -1);
00266 
00267    // calculate the p-value using the fast MCMC algorithm
00268    double pvalue;
00269    if (this -> CalculatePValueFast(this -> GetBestFitParameters(), pvalue))
00270       fPValue = pvalue;
00271    else
00272       BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::Fit() : Could not use the fast p-value evaluation.");
00273 
00274    // print summary to screen
00275    this -> PrintShortFitSummary();
00276 
00277    // no error
00278    return 1;
00279 }

int BCHistogramFitter::Fit (  )  [inline]

Performs the fit.

Returns:
An error code.

Definition at line 138 of file BCHistogramFitter.h.

00139          { return this -> Fit(fHistogram, fFitFunction); };

double BCHistogramFitter::FitFunction ( std::vector< double >  x,
std::vector< double >  parameters 
) [virtual]

Plots the histogram

Parameters:
options Options for plotting.
filename Name of the file which the histogram is printed into. The following options are available:
F : plots the fit function on top of the data E0 : plots the fit function and the 68% prob. uncertainty band of the fit function on top of the data E1 : plots the expectation from the fit function and the uncertainty bin-by-bin as error bars. Returns the y-value of the 1-dimensional fit function at an x and for a set of parameters.
x A vector with the x-value.
parameters A set of parameters.

Reimplemented from BCIntegrate.

Definition at line 230 of file BCHistogramFitter.cxx.

00231 {
00232    // set the parameters of the function
00233    fFitFunction -> SetParameters(&params[0]);
00234 
00235    return fFitFunction -> Eval(x[0]) * fHistogram -> GetBinWidth( fHistogram -> FindBin(x[0]) );
00236 }

TGraph* BCHistogramFitter::GetErrorBand (  )  [inline]
Returns:
pointer to the error band

Definition at line 72 of file BCHistogramFitter.h.

00073          { return fErrorBand; };

TF1* BCHistogramFitter::GetFitFunction (  )  [inline]
Returns:
The fit function

Definition at line 67 of file BCHistogramFitter.h.

00068          { return fFitFunction; };

TGraph* BCHistogramFitter::GetGraphFitFunction (  )  [inline]
Returns:
pointer to a graph for the fit function

Definition at line 77 of file BCHistogramFitter.h.

00078          { return fGraphFitFunction; };

TH1D* BCHistogramFitter::GetHistogram (  )  [inline]
Returns:
The histogram

Definition at line 62 of file BCHistogramFitter.h.

00063          { return fHistogram; };

double BCHistogramFitter::LogAPrioriProbability ( std::vector< double >  parameters  )  [virtual]

The log of the prior probability. Overloaded from BCModel.

Parameters:
parameters A vector of doubles containing the parameter values.

Reimplemented from BCModel.

Definition at line 152 of file BCHistogramFitter.cxx.

00153 {
00154    // using flat probability in all parameters
00155    double logprob = 0.;
00156    for(unsigned int i=0; i < this -> GetNParameters(); i++)
00157       logprob -= log(this -> GetParameter(i) -> GetRangeWidth());
00158 
00159    return logprob;
00160 }

double BCHistogramFitter::LogLikelihood ( std::vector< double >  parameters  )  [virtual]

The log of the conditional probability. Overloaded from BCModel.

Parameters:
parameters A vector of doubles containing the parameter values.

Reimplemented from BCModel.

Definition at line 164 of file BCHistogramFitter.cxx.

00165 {
00166    // initialize probability
00167    double loglikelihood = 0;
00168 
00169    // set the parameters of the function
00170    fFitFunction -> SetParameters(&params[0]);
00171 
00172    // get the number of bins
00173    int nbins = fHistogram -> GetNbinsX();
00174 
00175    // get bin width
00176    double dx = fHistogram -> GetBinWidth(1);
00177 
00178    // get function value of lower bin edge
00179    double fedgelow = fFitFunction -> Eval(fHistogram -> GetBinLowEdge(1));
00180    
00181    // loop over all bins
00182    for (int ibin = 1; ibin <= nbins; ++ibin)
00183    {
00184       // get upper bin edge
00185       double xedgehi = fHistogram -> GetBinLowEdge(ibin) + dx;
00186 
00187       // get function value at upper bin edge
00188       double fedgehi = fFitFunction -> Eval(xedgehi);
00189 
00190       // get the number of observed events
00191       double y = fHistogram -> GetBinContent(ibin);
00192 
00193       double yexp = 0.;
00194 
00195       // use ROOT's TH1D::Integral method
00196       if (fFlagIntegration)
00197          yexp = fFitFunction -> Integral(xedgehi-dx, xedgehi);
00198 
00199       // use linear interpolation
00200       else
00201       {
00202          yexp = fedgelow * dx + (fedgehi - fedgelow) * dx / 2.;
00203 
00204          // make the upper edge the lower edge for the next iteration
00205          fedgelow = fedgehi;
00206       }
00207 
00208       // get the value of the Poisson distribution
00209       loglikelihood += BCMath::LogPoisson(y, yexp);
00210    }
00211 
00212 // // get bin boundaries
00213 // double xmin = fHistogram -> GetBinLowEdge(ibin);
00214 // double xmax = fHistogram -> GetBinLowEdge(ibin+1);
00215 
00216 // // get the number of observed events
00217 // double y = fHistogram -> GetBinContent(ibin);
00218 
00219 // // get the number of expected events
00220 // double yexp = fFitFunction -> Integral(xmin, xmax);
00221 
00222 // // get the value of the Poisson distribution
00223 // loglikelihood += BCMath::LogPoisson(y, yexp);
00224 
00225    return loglikelihood;
00226 }

int BCHistogramFitter::SetFitFunction ( TF1 *  func  ) 
Parameters:
func The fit function @ return An error code (1:pass, 0:fail).

Definition at line 109 of file BCHistogramFitter.cxx.

00110 {
00111    // check if function exists
00112    if(!func)
00113    {
00114       BCLog::Out(BCLog::error,BCLog::error,"BCHistogramFitter::SetFitFunction() : TF1 not created.");
00115       return 0;
00116    }
00117 
00118    // set the function
00119    fFitFunction = func;
00120 
00121    // update the model name to contain the function name
00122    this -> SetName(TString::Format("HistogramFitter with %s",fFitFunction->GetName()));
00123 
00124    // reset parameters
00125    fParameterSet -> clear();
00126 
00127    // get the new number of parameters
00128    int n = func -> GetNpar();
00129 
00130    // add parameters
00131    for (int i = 0; i < n; ++i)
00132    {
00133       double xmin;
00134       double xmax;
00135 
00136       func -> GetParLimits(i, xmin, xmax);
00137 
00138       this -> AddParameter(func->GetParName(i), xmin, xmax);
00139    }
00140 
00141    // no error
00142    return 1;
00143 }

void BCHistogramFitter::SetFlagIntegration ( bool  flag  )  [inline]

Sets the flag for integration.
true: use ROOT's TH1D::Integrate()
false: use linear interpolation

Definition at line 101 of file BCHistogramFitter.h.

00102          { fFlagIntegration = flag; };

int BCHistogramFitter::SetHistogram ( TH1D *  hist  ) 
Parameters:
hist The histogram @ return An error code (1:pass, 0:fail).

Definition at line 55 of file BCHistogramFitter.cxx.

00056 {
00057    // check if histogram exists
00058    if (!hist)
00059    {
00060       BCLog::Out(BCLog::error,BCLog::error,"BCHistogramFitter::SetHistogram() : TH1D not created.");
00061       return 0;
00062    }
00063 
00064    // set pointer to histogram
00065    fHistogram = hist;
00066 
00067    // create a data set. this is necessary in order to calculate the
00068    // error band. the data set contains as many data points as there
00069    // are bins. for now, the data points are empty.
00070    BCDataSet * ds = new BCDataSet();
00071 
00072    // create data points and add them to the data set.
00073    int nbins = fHistogram -> GetNbinsX();
00074    for (int i = 0; i < nbins; ++i)
00075    {
00076       BCDataPoint* dp = new BCDataPoint(2);
00077       ds -> AddDataPoint(dp);
00078    }
00079 
00080    // set the new data set.
00081    this -> SetDataSet(ds);
00082 
00083    // calculate the lower and upper edge in x.
00084    double xmin = hist -> GetBinLowEdge(1);
00085    double xmax = hist -> GetBinLowEdge(nbins+1);
00086 
00087    // calculate the minimum and maximum range in y.
00088    double histymin = hist -> GetMinimum();
00089    double histymax = hist -> GetMaximum();
00090 
00091    // calculate the minimum and maximum of the function value based on
00092    // the minimum and maximum value in y.
00093    double ymin = TMath::Max(0., histymin - 5.*sqrt(histymin));
00094    double ymax = histymax + 5.*sqrt(histymax);
00095 
00096    // set the data boundaries for x and y values.
00097    this -> SetDataBoundaries(0, xmin, xmax);
00098    this -> SetDataBoundaries(1, ymin, ymax);
00099 
00100    // set the indeces for fitting.
00101    this -> SetFitFunctionIndices(0, 1);
00102 
00103    // no error
00104    return 1;
00105 }

Member Data Documentation

TGraph* BCHistogramFitter::fErrorBand [private]

Pointer to the error band (for legend)

Definition at line 179 of file BCHistogramFitter.h.

TF1* BCHistogramFitter::fFitFunction [private]

The fit function

Definition at line 170 of file BCHistogramFitter.h.

bool BCHistogramFitter::fFlagIntegration [private]

Flag for using the ROOT TH1D::Integral method (true), or linear interpolation (false)

Definition at line 175 of file BCHistogramFitter.h.

TGraph* BCHistogramFitter::fGraphFitFunction [private]

Pointer to a graph for displaying the fit function

Definition at line 183 of file BCHistogramFitter.h.

TH1D* BCHistogramFitter::fHistogram [private]

The histogram containing the data.

Definition at line 166 of file BCHistogramFitter.h.

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