In [1]:
import ROOT
ROOT.RooMsgService.instance().setGlobalKillBelow(5)
%jsroot on

from matplotlib import pyplot
%matplotlib inline

Welcome to JupyROOT 6.07/07

Reading Workspaces and Recreating Models

First we read in the signal region model from the workspace building notebook.

In [2]:
infile = ROOT.TFile.Open("output/SRselection.root")
w = infile.Get("w")
mc = w.obj("ModelConfig")
data = w.data("data")
mc.LoadSnapshot()
x = w.var("x")

Bayesian Markov Chain Monte Carlo Calculator

RooStats has several built in Bayesian calculators. We have the choice between using the MCMC calculator or the numerical one. We will start using the MCMC and optionally in the next exercise we will show how to use the numerical Bayesian calculator. We will use the same model used before for the Profile Likelihood calculator and repeat a similar process.

  • Read the model from the file as in the previous exercise
  • Create the MCMCCalculator class
  • Configure the calculator setting the proposal function, number of iterations, etc..
  • Call the GetInterval function to compute the desired interval. The returned MCMCInterval class will contained the resulting Markov-Chain.
  • Plot the result using the MCMCIntervalPlot class

Once the Data etc are extracted from the workspace. We also need the number of signal events.
To optimise we take a range of possible values (Truth was 100 events)

In [3]:
nsig = w.var("n_sig")
nsig.setRange(300,600)

Now we define our MCMC calculator. Again looking at 68.3% confidence level.

In [4]:
mcmc = ROOT.RooStats.MCMCCalculator(data,mc)
mcmc.SetConfidenceLevel(0.683)

Now MCMC methods have various implementations and approaches. Here we must propose a function. The sequential proposal given is usually fairly robust.

In [5]:
sp = ROOT.RooStats.SequentialProposal(0.1)
mcmc.SetProposalFunction(sp)

If you're interested in this I'd recommend looking it up. Here I just define lots of parameters. SetNumIters is related to the number of iterations of the Metropolis-Hastings algorithm SetBurnInSteps concerns how many interations are needed to avoid local minima due to low stats We can construct either a one-sided Bayesian interval (SetLetSideTailFraction=0) or central bayesian interval (SetLetSideTailFraction=0.5)

In [6]:
mcmc.SetNumIters(100000)
mcmc.SetNumBurnInSteps(1000)
mcmc.SetLeftSideTailFraction(0.5)

Now we run the calculator! (may take a while depending on the number of steps)

In [7]:
interval = mcmc.GetInterval()

Metropolis-Hastings progress: ....................................................................................................

Now we plot and print!

In [8]:
firstPOI = mc.GetParametersOfInterest().first()

print "68% interval on {0} is : [{1},{2}]".format(firstPOI.GetName(),interval.LowerLimit(firstPOI),interval.UpperLimit(firstPOI))

68% interval on n_sig is : [346.779036895,527.144738904]
In [9]:
c = ROOT.TCanvas()
plot = ROOT.RooStats.MCMCIntervalPlot(interval)
plot.Draw()
c.Draw()