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 0.9

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₁ = 50 is the number of background only counts, and N₂ = 55 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 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.9/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. The makefile is custom fit to the installation on your system, and in the following solutions, we will not display it.

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 "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, N_1)=(100, 50)
   BCDataPoint * backgroundMeasurement = new BCDataPoint(2);
   backgroundMeasurement->SetValue(0, 100);
   backgroundMeasurement->SetValue(1, 50);

   // 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");

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

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

   return 0;
}

In file CountingExp.cxx

#include <BAT/BCMath.h>

//...

void CountingExp::DefineParameters()
{
   // Allowed range for R_B is [0, 2]
   AddParameter("R_B", 0.0, 2.0);
}

double CountingExp::LogLikelihood(const 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(const std::vector<double> & parameters)
{
   // This method returns the logarithm of the prior probability for the
   // parameters p(parameters).

   double logprob = 0.;


   // normalize flat prior with parameter ranges
   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 0.5.

| Source code | Results | Plots |

Hide solution

Step 3 - Include the second measurement

Now add the second measurement which includes the source, and learn 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 <BAT/BCLog.h>
#include <BAT/BCAux.h>
#include <BAT/BCDataSet.h>
#include <BAT/BCH1D.h>

#include <TCanvas.h>

#include "CountingExp.h"

int main()
{

// ...


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

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

// ...

   // 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;

}

In file CountingExp.cxx

void CountingExp::DefineParameters()
{
   // Allowed range for R_B is [0, 2]
   AddParameter("R_B", 0.0, 5.0);

   // Allowed range for R_S is [0, 1]
   AddParameter("R_S", 0, 1.0);
}

double CountingExp::LogLikelihood(const 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. Verify that the variance of the posterior of R_B is smaller than in step 2, i.e. the additional measurement has constrained the value of R_B to a smaller region.

| Source code | Results | 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 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,1], R_S to the range [0,0.5] and use more bins (300 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 significantly less 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.0, 1.0);
   AddParameter("R_S", 0.0, 0.5);
}

double CountingExp::LogLikelihood(const 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(const 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.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(signalMeasurement);

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

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

   // perform your analysis here


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


   ...

   // 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());

   ...



   // prepare a new plot
   mycanvas->Clear();

   // 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");


   delete mycanvas;
   delete marginalHist1;
   delete legend;

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

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

   delete dataSet;
   delete m;

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

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

   return 0;
}

| Source code | Plot |

Hide solution

Frederik Beaujean

Last modified: Tue Apr 17 15:59:33 CEST 2012