[RF] RooStats.BayesianCalculator calculates the credibility interval on prior rather than on posterior

[capriotti]

Check duplicate issues.

• Checked for duplicates

Description

The RooStats::BayesianCalculator class seems to calculate the credibility interval based on the the prior PDF rather than on the posterior PDF unless RooStats::BayesianCalculator::SetScanOfPosterior(N) is explicitly called.

Reproducer

#import some modules
from ROOT import RooRealVar, RooFormulaVar, RooGaussian, RooExponential, RooAddPdf, RooFit
from ROOT import RooMinimizer, RooDataHist, RooDataSet, RooPlot, RooGenericPdf, RooArgSet, RooArgList
from RooFit import *

#declare the observable
mass = RooRealVar("mass","mass",6500,7500)

#make a toy MC: exponential
slope_gen = RooRealVar("slope_gen","slope_gen",-0.0008)
expo_gen = RooExponential("expo_gen","expo_gen",mass,slope_gen)

#declare the number of generated events
ngen = 1e3

#generate ngen events according to the total PDF
data = expo_gen.generate(RooArgSet(mass), ngen)

#declare the observable
mass = RooRealVar("mass","mass",6500,7500)

#create the model: Tbc signal
mean = RooRealVar("mean","mean",7024,6800,7200)
sigma = RooRealVar("sigma","sigma",44, 20, 60)
gauss = RooGaussian("gauss","gauss",mass,mean,sigma)

#create the model: background
slope = RooRealVar("slope","slope",-0.0008, -0.001, -0.0001)
expo = RooExponential("expo","expo",mass,slope)

#declare the signal and background yields
ngauss = RooRealVar("ngauss","ngauss",10,0,50)
nexpo = RooRealVar("nexpo","nexpo",1e3,5e2,2e6)

#sum the two PDFs to create the total gauss+expo PDF
pdftot = RooAddPdf("pdftot","pdftot",RooArgList(gauss,expo),RooArgList(ngauss,nexpo))

w = ROOT.RooWorkspace("ws")
w.Import(data)
w.Import(pdftot)
w.defineSet("poi", "ngauss")
w.defineSet("obs", "mass")
w.var("mean").setConstant(1)
w.var("sigma").setConstant(1)
w.var("slope").setConstant(1)
w.var("nexpo").setConstant(1)

# Create the ModelConfig for the signal + background hypothesis
sbModel = ROOT.RooStats.ModelConfig("sig+bkg model")
# Import the RooWorkspace
sbModel.SetWorkspace(w)
# Set the meaning of the workspace objects
sbModel.SetPdf(w["pdftot"])
sbModel.SetParametersOfInterest(w.set("poi"))
sbModel.SetObservables(w.set("obs"))

# Create the ModelConfig for the background only hypothesis
bModel = sbModel.Clone("bkg-only model")
# Set the meaning of the workspace objects
bModel.SetPdf(w["pdftot"])
poi = bModel.GetParametersOfInterest().first()
poi.setVal(0)
bModel.SetSnapshot(poi)

from ROOT import RooUniform
prior = RooUniform("prior","prior",ngauss)

#Construct the bayesian calculator
bc = ROOT.RooStats.BayesianCalculator(w.data("expo_genData"), sbModel)
#bc = ROOT.RooStats.BayesianCalculator(data, pdftot, {ngauss}, prior, 0)
bc.SetConfidenceLevel(0.95)
bc.SetIntegrationType("Gauss")
bc.SetLeftSideTailFraction(0.) # for upper limit
#bc.SetScanOfPosterior(50)
#
bcInterval = bc.GetInterval()

ROOT version

root_v6.32.04.Linux-ubuntu22.04-x86_64-gcc11.4

Installation method

pre-built binary

Operating system

Ubuntu 22.04.3

Additional context

This has been tested on Google Colab using python 3.11