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.

Charged current cross-section analysis - tutorial for BAT 0.9

Physics motivation

The charged current (CC) deep inelastic scattering cross-section in e⁺p interactions, σ_CC, was measured with a polarized positron beam provided by the HERA accelerator. The Standard Model predicts that the CC cross-section depends linearly on the polarization P, i.e., σ_CC = σ_CC(P) ∝ P, and vanishes for a polarization of P = -1 due to the (V-A) coupling of the W-boson to left-handed particles. A limit on a right-handed contribution to the cross-section is set in the following analysis, i.e., on the quantity σ_CC(-1). Further information on the details of the physics and the original analysis can be found here.

Tutorial

This tutorial shows how to implement a Bayesian model into BAT. In the first part a straight line fit is done in order to estimate the limit on the parameter under study. Subsequently, it is shown how prior knowledge from previous measurements can be used to improve the results. Sources of systematic uncertainties are taken into account as nuisance parameters in the third part. Finally, the result obtained is combined with a second meausrement of cross-section.

The tutorial is split into seven steps:

Step 1 - Getting started
Step 2 - Reading in the data
Step 3 - Defining the model
Step 4 - Calculating the limit
Step 5 - Including prior knowledge *
Step 6 - Including systematic uncertainties *
Step 7 - Combining results *

Steps marked with * are advanced examples and independent of the limit part of the tutorial. A proposal for additional studies is also given.

Step 1 - Getting started

Run the script CreateProject.sh in the tools directory to create a new project (CCXSecAnalysis) and a model class (CCXSec). Read through the code and compile it.

Solution

Step 2 - Reading in the data

The e⁺p data is stored in an ASCII file with three columns: polarization (P), measured cross-section (σ) and (symmetric) uncertainty on the measured cross-section (Δσ). The uncertainties on each measurement are assumed to be Gaussian. The file can be found here.

Modify the runCCXSecAnalysis.cxx code to read in the data set. Use the method BCDataSet::ReadDataFromFile (find out more about the method in the introduction to BAT).

Draw the data set (σ vs. P). The data can be accessed by the method BCDataSet::GetDataPoint(i)->GetValue(j).

What do you expect (from the SM) for the cross-section at a polarization of P=-1?

Solution

In file runCCXSecAnalysis.cxx

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

#include <CCXSec.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 CCXSec object
  CCXSec * m = new CCXSec();

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

   // create a new summary tool object
   BCSummaryTool * summary = new BCSummaryTool(m);

  // read data set from file 
  BCDataSet * mydataset = new BCDataSet();
  mydataset->ReadDataFromFile("data.txt", 3);
  m->SetDataSet(mydataset); 

   // perform your analysis here

   // normalize the posterior, i.e. integrate posterior
   // over the full parameter space
//   m->Normalize();

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

   // run mode finding; by default using Minuit
//   m->FindMode();

   // 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("CCXSec_plots.ps");

   // print all summary plots
//   summary->PrintParameterPlot("CCXSec_parameters.eps");
//   summary->PrintCorrelationPlot("CCXSec_correlation.eps");
//   summary->PrintKnowledgeUpdatePlots("CCXSec_update.ps");

   // calculate p-value
//   m->CalculatePValue( m->GetBestFitParameters() );

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


  delete m;
  delete summary;

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

  // return without error 
  return 0;
}

In file runCCXSecAnalysis.cxx

#include <BAT/BCLog.h>
#include <BAT/BCAux.h>
#include <BAT/BCSummaryTool.h>
#include <BAT/BCDataSet.h>
#include <TGraphErrors.h> 
#include <TH2D.h>
#include <TCanvas.h> 

#include <CCXSec.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 CCXSec object
  CCXSec * m = new CCXSec();

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

   // create a new summary tool object
   BCSummaryTool * summary = new BCSummaryTool(m);

  // read data set from file 
  BCDataSet * mydataset = new BCDataSet();
  mydataset->ReadDataFromFile("data.txt", 3);
  m->SetDataSet(mydataset);

  // define histogram for data axes 
  TH2D * myhistaxes = new TH2D("myhistaxes", ";Polarization;#sigma_{CC} [pb]", 1, -1., 1., 1, 0.0, 80.); 
  myhistaxes->SetStats(kFALSE); 

  // copy data into graph
  TGraphErrors * mygraphdata = new TGraphErrors(mydataset->GetNDataPoints());
  mygraphdata->SetMarkerStyle(20);
  for (int i = 0; i < int(mydataset->GetNDataPoints()); ++i)
  {
    double P = mydataset->GetDataPoint(i)->GetValue(0); 
    double sig = mydataset->GetDataPoint(i)->GetValue(1); 
    double errsig = mydataset->GetDataPoint(i)->GetValue(2); 

    mygraphdata->SetPoint(i, P, sig);
    mygraphdata->SetPointError(i, 0.0, errsig);
  } 
  
  // print data 
  TCanvas * mycanvas = new TCanvas("mycanvas"); 
  mycanvas->cd();
  myhistaxes->Draw(); 
  mygraphdata->Draw("SAMEPE"); 
  mycanvas->Print("data.eps"); 

   // perform your analysis here

   // normalize the posterior, i.e. integrate posterior
   // over the full parameter space
//   m->Normalize();

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

   // run mode finding; by default using Minuit
//   m->FindMode();

   // 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("CCXSec_plots.ps");

   // print all summary plots
//   summary->PrintParameterPlot("CCXSec_parameters.eps");
//   summary->PrintCorrelationPlot("CCXSec_correlation.eps");
//   summary->PrintKnowledgeUpdatePlots("CCXSec_update.ps");

   // calculate p-value
//   m->CalculatePValue( m->GetBestFitParameters() );

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


  delete m;
  delete summary;

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

  // return without error 
  return 0;
}

The plot of the data can be found here.

The cross-section for a polarized positron beam is expected to be 0 due to the V-A structure of the weak interaction.

| Source code |

Hide solution

Step 3 - Defining the model

It is experimentally difficult to reach a polarization of more than 50%. However, the Standard Model predicts a linear dependence of the cross-section on the polarization. This dependence can be used to extrapolate the cross-section to a polarization of P=-1. In the following, the data is fitted with a straight line. The two parameters of the line are chosen to be the cross-section at a polarization of P=0, σ₀, and the cross-section at a polarization of P=-1, σ₁. The cross-section can be parameterized as

σ_CC(P) = P (σ₀ - σ₁) + σ₀

Using the Bayesian approach, the posterior probability density for the two parameters given the data set,

p(σ₀, σ₁ | data) ∝ p(data|σ₀, σ₁) ⋅ p₀(σ₀, σ₁)

is maximized. The measurements at different polarizations are assumed to be independent. The likelihood can then be a product of Gaussian terms. Modify the model code accordingly:

Add the two parameters to the model. No prior knowledge about the parameters is used in the fit, i.e., p₀(σ_0/1) = const. Choose the ranges of the parameters to be larger than zero. Why?

Modify the method CCXSec::LogLikelihood(). The Gaussian terms are added logarithmically (mind the sign!). One can use the method BCMath::LogGaus()

Run the Markov Chain. Call the method BCModel::MarginalizeAll() in the file runCCXSecAnalysis.cxx.

Create the marginalized distributions of the parameters. Call the method BCModel::PrintAllMarginalized("plots.ps") in the file runCCXSecAnalysis.cxx. This will create all possible 1-D and 2-D marginalized distributions and write them in the file plots.ps.

Print the marginalized distribution including the marked 95% probability region using BCModel::GetMarginalized("sig1")->Draw(0, -95).

Write out the results into a .txt file. Use the method BCModel::PrintResults(const char * file).

Solution

In file CCXSec.cxx

void CCXSec::DefineParameters()
{
  // Add parameters to your model here.
  // You can then use them in the methods below by calling the
  // parameters.at(i) or parameters[i], where i is the index
  // of the parameter. The indices increase from 0 according to the
  // order of adding the parameters.

  AddParameter("sig0", 0.0, 50.0); 
  AddParameter("sig1", 0.0, 20.0); 
}

Cross-section is a physical quantity which has to be larger or equal to zero.

In file CCXSec.cxx


double CCXSec::LogLikelihood(const std::vector  & 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 parameters 
  double sig0 = parameters.at(0); 
  double sig1 = parameters.at(1); 

  for (int i = 0; i < GetNDataPoints(); ++i)
  {
    // get data point 
    double P      = GetDataPoint(i)->GetValue(0); 
    double sig    = GetDataPoint(i)->GetValue(1); 
    double errsig = GetDataPoint(i)->GetValue(2); 

    // calculate expected cross-section 
    double sigexp = P * (sig0 - sig1) + sig0; 

    // calculate likelihood 
    logprob += BCMath::LogGaus(sig, sigexp, errsig);
  }

  return logprob;   
}

In file runCCXSecAnalysis.cxx

int main()
{
  ...
  // run MCMC 
  m->MarginalizeAll();
  ...
}

In file runCCXSecAnalysis.cxx

int main()
{
  ...
  // run MCMC 
  m->MarginalizeAll();

  // print all marginalized distributions 
  m->PrintAllMarginalized("plots.ps"); 
  ...
}

In file runCCXSecAnalysis.cxx

#include <BAT/BCH1D.h> 

int main()
{
  ... 
  // run MCMC 
  m->MarginalizeAll();

  // print all marginalized distributions 
  m->PrintAllMarginalized("plots.ps"); 

  // print marginalized distribution for sig1 including 95% limit
  mycanvas->cd();
  m->GetMarginalized("sig1")->Draw(0, -95);
  mycanvas->Print("sig1.eps");
  ...
}

In file runCCXSecAnalysis.cxx

#include <BAT/BCH1D.h> 

int main()
{
  ... 
  // run MCMC 
  m->MarginalizeAll();

  // print all marginalized distributions 
  m->PrintAllMarginalized("plots.ps"); 

  // print marginalized distribution for sig1 including 95% limit
  mycanvas->cd();
  m->GetMarginalized("sig1")->Draw(0, -95);
  mycanvas->Print("sig1.eps"); 

  // print results 
  m->PrintResults("results.txt");
  ...
}

| Source code |

Hide solution

Step 4 - Calculating the limit

Compile and run the code. Look at the 1-d and 2-d distributions. What is the 95% limit on the cross-section at a polarization of -1?

Note that the Markov Chains run for 100,000 iterations by default. Change this number in order to decrease the statistical uncertainty on the limit. Use the method BCModel::MCMCSetNIterationsRun(int n).

Solution

Step 5 - Including prior knowledge

The CC cross-section has previously been measured with a beam of unpolarized positrons. The measurement yielded σ₀ = (36.2^+1.1_-0.49) pb. This result can be used as a prior for the parameter σ₀ in the current analysis.

Include the prior knowledge on the parameter σ₀ into the model. The probability density can be parameterized by the following function
p₀(σ₀)=(σ₀- 35.0 pb)⋅e^{-(σ₀- 35.0} pb)²/2.88 pb² with σ₀ > 35.0 pb.

Repeat step number 4 and compare the limits obtained.

Solution

Step 6 - Including systematic uncertainties

So far, the uncertainties on the model parameters are calculated by propagating the uncertainties from the individual measurements. Additional systematic uncertainties can be included as nuisance parameters which have to be added to the model:

The parameter δ_p represents the uncertainty on the polarization measurement. A Gaussian distribution centered around 0 with a width of 0.02 is used to parameterize the uncertainty. P is then replaced by P(1+δ_p) in the likelihood. The Gaussian is used as a prior for δ_p.

The parameter δ_s represents all other sources of systematic uncertainty. A flat distribution centered around 0 with a width of 0.034 is used to parameterize the uncertainty. The cross-section is then multiplied by (1+δ_s) in the likelihood.

Include the extra nuisance parameters dp and ds.

Add the prior knowledge about the parameters to the model.

Repeat step number 4 and compare the limit obtained.

Solution

In file CCXSec.cxx

void CCXSec::DefineParameters()
{
  // Add parameters to your model here.
  // You can then use them in the methods below by calling the
  // parameters.at(i) or parameters[i], where i is the index
  // of the parameter. The indices increase from 0 according to the
  // order of adding the parameters.

  AddParameter("sig0",   0.0,    50.0); 
  AddParameter("sig1",   0.0,    20.0); 
  AddParameter("dp",  -  0.1,     0.1);
  AddParameter("ds",  -  0.017,   0.017);
}

In file CCXSec.cxx

double CCXSec::LogAPrioriProbability(const std::vector  & parameters)
{
  // This method returns the logarithm of the prior probability for the
  // parameters p(parameters).

  double logprob = 0.;

  // get parameters 
  double sig0 = parameters.at(0);
  double dp   = parameters.at(2);

  // use result from measurement with unpolarized beam as prior for
  // parameter sig0
  if (sig0 > 35.)
  logprob += log(sig0 - 35.0) - (sig0 - 35.0) * (sig0 - 35.0) / 2.88;
  else
  logprob += -999;  

  // use prior knowledge of nuisance parameter dp
  logprob += BCMath::LogGaus(dp, 0.0, 0.02);

  return logprob;
}

and

double CCXSec::LogLikelihood(const std::vector  & 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 parameters 
  double sig0 = parameters.at(0); 
  double sig1 = parameters.at(1); 
  double dp   = parameters.at(2);
  double ds   = parameters.at(3); 

  for (int i = 0; i < GetNDataPoints(); ++i)
  {
    // get data point 
    double P      = GetDataPoint(i)->GetValue(0); 
    double sig    = GetDataPoint(i)->GetValue(1); 
    double errsig = GetDataPoint(i)->GetValue(2); 

    // calculate expected cross-section 
    double sigexp = (P * (1 + dp) * (sig0 - sig1) + sig0) * (1 + ds); 

    // calculate likelihood 
    logprob += BCMath::LogGaus(sig, sigexp, errsig);
  }

  return logprob;   
}

The limit on the parameter σ₁ is 5.66 pb. The full posterior probability density for the parameter can be found here.

| Source code | Results | 1D and 2D Marginalized plots |

Hide solution

Step 7 - Combining results

A second experiment measures the cross-section at a polarization of P=-1. The posterior for the parameter σ₁ from the second experiment has the shape of a Gaussian centered around 0.5 pb with a width of 4 pb, and is zero for values smaller than 0 pb. Combine both results and calculate the 95% limit on σ₁.

Solution

Additional studies

Try different priors on the model parameters. E.g., a theory could predict an expotential for σ₁. How does the limit on the parameter change?

Try different priors on the nuisance parameters. How does the limit change if the prior for δ_p is assumed to be flat?

Calculate the uncertainty band on the fit function.

Study the effect on the limit on σ₁ if it is not constrained to positive values.

Further material and documentation

The original analysis can be found here:

J. Sutiak,
Measurement of the Polarised Charged Current Cross Sections with the ZEUS Detector,
PhD thesis, TU Munich, 2006.

Current results from ZEUS and H1:

S. Chekanov et al.,
Measurement of Charged Current Deep Inelastic Scattering Cross Sections with a Longitudinally Polarised Electron Beam at HERA,
Eur. Phys. J. C 61 (2009) 223

A. Aktas et al.,
First Measurement of Charged Current Cross Sections at HERA with Longitudinally Polarised Positrons,
Phys. Lett. B 634 (2006) 173 [hep-ex/0512060]

For an overview on the V-A structure of the weak interaction, see e.g.,

F. Halzen, A. Martin,
Quarks & Leptons,
John Wiley & Sons, 1984, Ch. 12

D. Griffiths,
Introduction to Elementary Particles,
Wiley-VCH, 2nd edition, 2008, Ch. 9

For a short introduction to this tutorial, see this PDF.

Kevin Kröninger

Last modified: Thu Oct 6 03:27:34 CEST 2011