BAT - Bayesian Analysis Toolkit : tutorials

This C++ version of BAT is still being maintained, but addition of new features is unlikely. Check out our new incarnation, BAT.jl, the Bayesian analysis toolkit in Julia. In addition to Metropolis-Hastings sampling, BAT.jl supports Hamiltonian Monte Carlo (HMC) with automatic differentiation, automatic prior-based parameter space transformations, and much more. See the BAT.jl documentation.

Measuring a decay rate - tutorial for BAT

Physics motivation

When developing a new detector in particle or nuclear physics it is important to measure a detector's efficiency. In practice this is often accomplished in two steps. First the background of the environment is measured. Then a radiation source of known strength is placed near the detector. Finally the counts with and without source are compared to verify that the detector operates properly.

In the following we assume that two measurements have been made, both of duration T = 100s. N₁ = 100 is the number of background only counts, and N₂ = 110 is the number of counts including the source. N₁ and N₂ are supposed to arise from a Poisson process. The goal is to extract the background rate R_B and the signal rate R_S.

An expanded introduction including the necessary formulae for this tutorial is found /Introduction_Bologna.pdf">here.

Tutorial

This tutorial introduces the basics of BAT. It shows how to

create a simple model
compile and run the code
extract best-fit parameters and associated uncertainties
plot marginalized distributions and parameter correlations

We assume the reader has downloaded and installed BAT and is familiar with the C++ programming language. The detailed reference guide documenting BAT's internals provides further assistance.

The tutorial is split into four steps:

Step 1 - Compiling your first BAT program
Step 2 - Fitting the background-only model
Step 3 - Including the signal contribution
Step 4 - The next steps

Step 1 - Getting started

Make sure the installation of BAT was successful:

Navigate to the tools subdirectory of your BAT installation, e.g. ~/BAT-0.4.2/tools and run the script CreateProject.sh to create a project named CountingExp. A subdirectory is created.
Have a look at the generated C++ classes and compile the code using make.

Solution

Step 2 - Fitting the background-only model

The first model we test considers only background.

Create a data point, add it to a data set and register the data set with the model
Define the parameter R_B and add it to the model
Define the log likelihood for the Poisson process with parameter R_B. The natural logarithm of the factorial is provided by BCMath::LogFact(int n). One can also use the approximation provided by BCMath::ApproxLogFact(int n) which is much faster for large numbers.
Use a flat prior for R_B
Start to sample from the posterior using the Markov chain
Find the mode of the posterior
Save the results of the fit in text form and create a plot of the (marginal) posterior distribution

Solution

To execute your code, type:

$> make $> ./runCountingExp

In file runCountingExp.cxx

#include <BAT/BCLog.h>
#include <BAT/BCAux.h>
#include <BAT/BCDataSet.h>
#include <BAT/BCDataPoint.h>

#include "CountingExp.h"

int main()
{
   // set nice style for drawing than the ROOT default
   BCAux::SetStyle();

   // open log file
   // BCLog::OpenLog("log.txt");
   BCLog::SetLogLevel(BCLog::detail);

   // create new CountingExp object
   CountingExp * m = new CountingExp();
   //create background measurement data point (T,R_B)=(100, 100)
   BCDataPoint * backgroundMeasurement = new BCDataPoint(2);
   backgroundMeasurement->SetValue(0, 100);
   backgroundMeasurement->SetValue(1, 100);

   // add the single single measurement to the data set
   BCDataSet * dataSet = new BCDataSet();
   dataSet->AddDataPoint(backgroundMeasurement);

   // register the data set with the model
   m->SetDataSet(dataSet);

   BCLog::OutSummary("Test model created");

   // perform your analysis here


   // run MCMC and marginalize posterior wrt. all parameters
   // and all combinations of two parameters
   m->MarginalizeAll();

   // if MCMC was run before (MarginalizeAll()) it is
   // possible to use the mode found by MCMC as
   // starting point of Minuit minimization
   m->FindMode(m->GetBestFitParameters());

   // draw all marginalized distributions into a PostScript file
   m->PrintAllMarginalized("CountingExp_plots.ps");

   // print results of the analysis into a text file
   m->PrintResults("CountingExp_results.txt");

   // close log file
   BCLog::CloseLog();


   delete backgroundMeasurement;
   delete dataSet;
   delete m;

   BCLog::OutSummary("Test program ran successfully");
   BCLog::OutSummary("Exiting");

   return 0;
}

In file CountingExp.cxx

#include <BAT/BCMath.h>

//...

void CountingExp::DefineParameters()
{

   // Allowed range for R_B is [0, 5]
   AddParameter("R_B", 0.0, 5.0);
}

double CountingExp::LogLikelihood(std::vector <double> parameters)
{
   // This methods returns the logarithm of the conditional probability
   // p(data|parameters). This is where you have to define your model.

   double logprob = 0.;


   // get background measurement
   double T = GetDataPoint(0)->GetValue(0);
   double N1 = GetDataPoint(0)->GetValue(1);

   // extract value of background rate
   double R_B = parameters.at(0);

   // calculate expected counts given background rate
   double n_B = R_B * T;

   // update likelihood
   logprob += -n_B + N1 * log(n_B) - BCMath::LogFact(N1);

   return logprob;
}

double CountingExp::LogAPrioriProbability(std::vector<double> parameters)
{
   // This method returns the logarithm of the prior probability for the
   // parameters p(parameters).

   double logprob = 0.;


   // For flat prior it's very easy.
   for (unsigned int i = 0; i < GetNParameters(); i++)
      logprob -= log(GetParameter(i)->GetRangeWidth());

   return logprob;
}

Hint: the best estimate for R_B should be around 1.0 . The posterior is not normalized in this example to save computing time, the value of mode, mean... of R_B are unaffected by the normalization.

| /step2.tgz">Source code | /step2/CountingExp_results.txt">Results | /step2/CountingExp_plots.ps">Plots |

Hide solution

Step 3 - Include the second measurement

Now analyze what can be learned from the second measurement with the source about signal and background rate.

Add a second data point, N₂, to the data set
Include the second parameter R_S with flat prior in the model
Update the likelihood to incorporate N₂ and R_S
Plot the marginal distributions and compare the values of mean, median and mode for the individual parameters. What is the correlation between R_B and R_S?
Extract the 95% limit on R_S and save the plot.

Solution

In file runCountingExp.cxx

#include <TCanvas.h>
#include <BAT/BCH1D.h>
#include "CountingExp.h"


int main()
{

// ...


   // create background + signal data point (T,N_2)=(100, 110)
   BCDataPoint * bsMeasurement = new BCDataPoint(2);
   bsMeasurement->SetValue(0, 100);
   bsMeasurement->SetValue(1, 110);

   // add the single measurement to the data set
   BCDataSet * dataSet = new BCDataSet();
   dataSet->AddDataPoint(backgroundMeasurement);
   dataSet->AddDataPoint(bsMeasurement);

// ...

   // draw all marginalized distributions into a PostScript file
   m->PrintAllMarginalized("CountingExp_plots.ps");

   TCanvas * mycanvas = new TCanvas("mycanvas");
   m->GetMarginalized("R_S")->Draw(0, -95);
   mycanvas->Print("Limit.eps");
   delete mycanvas;

//..

   delete backgroundMeasurement;
   delete bsMeasurement;
   delete dataSet;
   delete m;

}

In file CountingExp.cxx

void CountingExp::DefineParameters()
{
   // allowed ranges [0, 5] for R_B and R_S
   AddParameter("R_B", 0.0, 5.0); // index 0

   AddParameter("R_S", 0, 5.0);   // index 1
}

double CountingExp::LogLikelihood(std::vector<double> parameters)
{
   // This methods returns the logarithm of the conditional probability
   // p(data|parameters). This is where you have to define your model.

   double logprob = 0.;


   // get background measurement
   double T1 = GetDataPoint(0)->GetValue(0);
   double N1 = GetDataPoint(0)->GetValue(1);

   // get background + signal measurement
   double T2 = GetDataPoint(1)->GetValue(0);
   double N2 = GetDataPoint(1)->GetValue(1);

   // extract value of background rate
   double R_B = parameters.at(0);


   // extract value of background + signal rate
   double R_S = parameters.at(1);

   // calculate expected counts given background rate
   double n_B = R_B * T1;

   // calculate expected counts for second measurement
   double n = (R_B + R_S) * T2;

   // update likelihood for both measurements
   logprob += -n_B + N1 * log(n_B) - BCMath::LogFact(N1);
   logprob += -n + N2 * log(n) - BCMath::LogFact(N2);

   return logprob;
}

The global mode should be around R_B = 1, R_S = 0.1. Make sure the variance of the posterior of R_B is smaller than in step 2. The additional measurement has constrained the value of R_B to a smaller region.

| /step3.tgz">Source code | /step3/CountingExp_results.txt">Results | /step3/CountingExp_plots.ps">Plots |

Hide solution

Step 4 - Further steps

The basic steps should be clear by now. We want to learn about more detailed features of BAT. A frequent task is to create plots. But what how to customize the plots? BAT resorts to the plotting engine provided by ROOT. As a simple example we create a plot that combines the posterior distributions P(R_B|N₁) and P(R_B|N₁,N₂) from step 2 and 3. Another issue is speed: during the Markov chain sampling, the BCModel methods LogLikelihood and LogAPrioriProbability are called in each iteration. Accordingly, even a small speedup in these methods may lead to a huge speedup of the overall execution time.

Redo step 2 and 3. Save P(R_B|N₁) and P(R_B|N₁,N₂) as a ROOT TH1D histogram. Limit R_B to the range [0,2] and use more bins (500 instead of the default 100) to store the marginalized distribution.
Normalize and plot the two histograms.
Measure the time it takes to run the program.
Modify LogAPrioriProbability to do nothing else than returning zero. This amounts to setting the prior to 1. Compare execution time.
Multiplying the likelihood by a constant just affects the normalization, but not the values of mode, mean... Thus remove all terms that are added to LogLikelihood and which are independent of R_B, R_S. You should observe that running the program takes only about a quarter of the time compared to 3.
Redo step 2, but now use the Reference prior (=Jeffrey's prior here) for R_B which reads P(R_B)∝1/√R_B. Does the posterior P(R_B|N₁) change significantly?

Solution

A typical shell session might give:

$#step 3.
$ time ./runCountingExp
real	0m15.042s
user	0m14.865s
sys	0m0.068s

$#step 4.
$ time ./runCountingExp
real	0m14.263s
user	0m14.213s
sys	0m0.048s

$#step 5.
$ time ./runCountingExp
real	0m2.458s
user	0m2.388s
sys	0m0.052s

In file CountingExp.cxx

void CountingExp::DefineParameters()
{

   // allowed ranges [0, 2] for R_B and R_S
   AddParameter("R_B", 0., 2);
   AddParameter("R_S", 0., 2);
}

double CountingExp::LogLikelihood(std::vector<double> parameters)
{
   // This methods returns the logarithm of the conditional probability
   // p(data|parameters). This is where you have to define your model.

   double logprob = 0.;

   // get background measurement
   double T1 = GetDataPoint(0)->GetValue(0);
   double N1 = GetDataPoint(0)->GetValue(1);

   // get background + signal measurement
   double T2=0;
   double N2=0;
   if (GetNDataPoints() > 1) {
      T2 = GetDataPoint(1)->GetValue(0);
      N2 = GetDataPoint(1)->GetValue(1);
   }


   // extract value of background rate
   double R_B = parameters.at(0);

   // extract value of background + signal rate
   double R_S = parameters.at(1);

   // calculate expected counts given background rate
   double n_B = R_B * T1;


   // calculate expected counts for second measurement
   double n=0;
   if (GetNDataPoints() > 1)
      n = (R_B + R_S) * T2;

   // update likelihood for both measurements
   // additive constants irrelevant for finding the mode, etc.
   // but may slow down evaluation

   // slow
//   logprob += -n_B + N1 * log(n_B) - BCMath::LogFact(N1);
//   if (GetNDataPoints() > 1)
//      logprob += -n + N2 * log(n) - BCMath::LogFact(N2);

   // fast
   logprob += -n_B + N1 * log(n_B);
   if (GetNDataPoints() > 1)
      logprob += -n + N2 * log(n);


    return logprob;
}

double CountingExp::LogAPrioriProbability(std::vector<double> parameters)
{
   // This method returns the logarithm of the prior probability for the
   // parameters p(parameters).


   // Constant prior is a constant prior. The actual value only required
   // for normalization, but not for mode finding.

   // slow

//   double logprob = 0.;

//   for (unsigned int i = 0; i < GetNParameters(); i++)
//      logprob -= log(GetParameter(i)->GetRangeWidth());

//   return logprob;

   // fast
   return 0.;

}

In file runCountingExp.cxx

#include <TCanvas.h>

#include <TH1D.h>
#include <TLegend.h>

#include "CountingExp.h"

int main()
{

//...

   // add the single measurement to the data set
   BCDataSet * dataSet = new BCDataSet();
   dataSet->AddDataPoint(backgroundMeasurement);

//   dataSet->AddDataPoint(bsMeasurement);

   // register the data set with the model
   m->SetDataSet(dataSet);

   BCLog::OutSummary("Test model created");

   // perform your analysis here


   m->SetNbins("R_B", 500);
   m->SetNbins("R_S", 500);


   // run MCMC and marginalize posterior wrt. all parameters
   // and all combinations of two parameters
   m->MarginalizeAll();

   // if MCMC was run before (MarginalizeAll()) it is
   // possible to use the mode found by MCMC as
   // starting point of Minuit minimization
   m->FindMode(m->GetBestFitParameters());


   // get marginal distribution P(R_B | N1 )
   TH1D * marginalHist1 = m->GetMarginalized("R_B")->GetHistogram();

   // create a copy of the histogram
   marginalHist1 = (TH1D*) marginalHist1->Clone("mHist1");

   // now fit again with N1 and N2
   dataSet->AddDataPoint(bsMeasurement);
   m->MarginalizeAll();
   m->FindMode(m->GetBestFitParameters());

   // get marginal distribution P(R_B|{N1,N2})
   TH1D * marginalHist2 = m->GetMarginalized("R_B")->GetHistogram();

   // create a canvas to hold the comparison plot
   TCanvas * mycanvas = new TCanvas("mycanvas");

   // draw second histogram in red
   marginalHist2->SetLineColor(2);

   // superimpose both distributions, normalize for easier comparison
   marginalHist2->DrawNormalized();
   marginalHist1->DrawNormalized("SAME");

   // add a legend to the plot
   TLegend * legend = new TLegend(0.7, 0.7, 0.85, 0.85);
   legend->AddEntry(marginalHist1, "N1", "l");
   legend->AddEntry(marginalHist2, "N1 and N2", "l");
   legend->Draw();

   // save a copy of the plot
   mycanvas->Print("R_B_update.eps");

   // print results of the analysis into a text file
   m->PrintResults("CountingExp_results.txt");

   // close log file
   BCLog::CloseLog();

   // clean up the heap
   delete backgroundMeasurement;
   delete bsMeasurement;
   delete dataSet;
   delete m;

   delete mycanvas;
   delete marginalHist1;
   delete legend;

   BCLog::OutSummary("Test program ran successfully");
   BCLog::OutSummary("Exiting");

   return 0;
}

| /step4.tgz">Source code | /step4/R_B_update.eps">Plot |

Hide solution

Frederik Beaujean

Last modified: Tue Feb 22 12:19:12 CET 2011