This example is a generalization of the on/off problem. It's a common setup for SUSY searches. Imagine that one has two variables "x" and "y" (eg. missing ET and SumET), see figure. The signal region has high values of both of these variables (top right). One can see low values of "x" or "y" acting as side-bands. If we just used "y" as a sideband, we would have the on/off problem.
If tau is known, this model is sufficient, but often tau is not known exactly. So one can use low values of "x" as an additional constraint for tau. Note that this technique critically depends on the notion that the joint distribution for "x" and "y" can be factorized. Generally, these regions have many events, so it the ratio can be measured very precisely there. So we extend the model to describe the left two boxes... denoted with "bar".
One can further expand the model to account for the systematic associated to assuming the distribution of "x" and "y" factorizes (eg. that tau is the same for off/on and offbar/onbar). This can be done in several ways, but here we introduce an additional parameter rho, which so that one set of models will use tau and the other tau*rho. The choice is arbitrary, but it has consequences on the numerical stability of the algorithms. The "bar" measurements typically have more events (& smaller relative errors). If we choose
the product tau*rho will be known very precisely (~1/sqrt(bbar)) and the contour in those parameters will be narrow and have a non-trivial tau~1/rho shape. However, if we choose to put rho on the non/noff measurements (where the product will have an error ~1/sqrt(b))
, the contours will be more amenable to numerical techniques. Thus, here we choose to define:
Left in this way, the problem is under-constrained. However, one may have some auxiliary measurement (usually based on Monte Carlo) to constrain rho. Let us call this auxiliary measurement that gives the nominal value of rho "rhonom". Thus, there is a 'constraint' term in the full model: P(rhonom | rho). In this case, we consider a Gaussian constraint with standard deviation sigma.
Note, the covariance matrix of the parameters has large off-diagonal terms. Clearly s,b are anti-correlated. Similarly, since noffbar >> nonbar, one would expect bbar,tau to be anti-correlated.
This can be seen below.
Similarly, since tau*rho appears as a product, we expect rho,tau to be anti-correlated. When the error on rho is significantly larger than 1/sqrt(bbar), tau is essentially known and the correlation is minimal (tau mainly cares about bbar, and rho about b,s). In the alternate parametrization (bbar* tau * rho) the correlation coefficient for rho,tau is large (and negative).
␛[1mRooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby␛[0m
Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
All rights reserved, please read http://roofit.sourceforge.net/license.txt
[#0] WARNING:InputArguments -- The parameter 'sigma' with range [-1e+30, 1e+30] of the RooGaussian 'mcCons' exceeds the safe range of (0, inf). Advise to limit its range.
[#1] INFO:ObjectHandling -- RooWorkspace::import(wspace) importing dataset modelData
[#1] INFO:Minization -- createNLL: caching constraint set under name CONSTR_OF_PDF_model_FOR_OBS_noff:noffbar:non:nonbar:rhonom with 0 entries
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoGLobalFit - find MLE
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoMinimizeNLL - using Minuit / Migrad with strategy 1
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:Minization -- The following expressions will be evaluated in cache-and-track mode: (on,off,onbar,offbar,mcCons)
[#1] INFO:Minization --
RooFitResult: minimized FCN value: 16.2872, estimated distance to minimum: 1.21273e-07
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
b 8.3602e+01 +/- 1.39e+01
bbar 9.9301e+02 +/- 3.15e+01
rho 1.2783e+00 +/- 1.99e-01
s 5.5397e+01 +/- 1.78e+01
tau 4.9405e+00 +/- 1.72e-01
Bayesian Calc. only supports on parameter of interest
[#1] INFO:Minization -- createNLL picked up cached constraints from workspace with 0 entries
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:Minization -- The following expressions will be evaluated in cache-and-track mode: (on,off,onbar,offbar,mcCons)
**********
** 1 **SET PRINT 1
**********
**********
** 2 **SET NOGRAD
**********
PARAMETER DEFINITIONS:
NO. NAME VALUE STEP SIZE LIMITS
1 b 1.17958e+02 1.38657e+01 0.00000e+00 3.00000e+02
2 bbar 1.00111e+03 3.15031e+01 5.00000e+02 2.00000e+03
3 rho 9.28979e-01 1.98664e-01 0.00000e+00 2.00000e+00
4 s 1.20959e+01 1.78108e+01 0.00000e+00 1.00000e+02
MINUIT WARNING IN PARAMETR
============== VARIABLE4 BROUGHT BACK INSIDE LIMITS.
5 tau 4.89226e+00 1.71715e-01 3.00000e+00 7.00000e+00
**********
** 3 **SET ERR 0.5
**********
**********
** 4 **SET PRINT 1
**********
**********
** 5 **SET STR 1
**********
NOW USING STRATEGY 1: TRY TO BALANCE SPEED AGAINST RELIABILITY
**********
** 6 **MIGRAD 2500 1
**********
FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.
START MIGRAD MINIMIZATION. STRATEGY 1. CONVERGENCE WHEN EDM .LT. 1.00e-03
FCN=18.2144 FROM MIGRAD STATUS=INITIATE 20 CALLS 21 TOTAL
EDM= unknown STRATEGY= 1 NO ERROR MATRIX
EXT PARAMETER CURRENT GUESS STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 b 1.17958e+02 1.38657e+01 9.47845e-02 -1.99668e-02
2 bbar 1.00111e+03 3.15031e+01 4.45476e-02 -1.41126e-01
3 rho 9.28979e-01 1.98664e-01 2.00529e-01 -1.44935e-02
4 s 1.20959e+01 1.78108e+01 5.77645e-01 -2.24297e+00
5 tau 4.89226e+00 1.71715e-01 8.60896e-02 -8.62115e-02
ERR DEF= 0.5
MIGRAD MINIMIZATION HAS CONVERGED.
MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
COVARIANCE MATRIX CALCULATED SUCCESSFULLY
FCN=16.2872 FROM MIGRAD STATUS=CONVERGED 190 CALLS 191 TOTAL
EDM=3.60951e-06 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 b 8.35905e+01 1.38673e+01 7.30675e-05 2.40056e-02
2 bbar 9.93029e+02 3.15018e+01 5.19045e-05 -3.02691e-02
3 rho 1.27859e+00 1.98755e-01 1.58308e-04 1.62269e-02
4 s 5.54105e+01 1.78121e+01 6.70504e-04 3.43680e-04
5 tau 4.94037e+00 1.71699e-01 9.49080e-05 -2.33172e-02
ERR DEF= 0.5
EXTERNAL ERROR MATRIX. NDIM= 25 NPAR= 5 ERR DEF=0.5
1.930e+02 8.358e+01 -2.620e+00 -1.929e+02 -5.000e-01
8.358e+01 9.930e+02 1.248e-04 -8.355e+01 -4.940e+00
-2.620e+00 1.248e-04 4.008e-02 2.619e+00 -7.467e-07
-1.929e+02 -8.355e+01 2.619e+00 3.319e+02 4.998e-01
-5.000e-01 -4.940e+00 -7.467e-07 4.998e-01 2.955e-02
PARAMETER CORRELATION COEFFICIENTS
NO. GLOBAL 1 2 3 4 5
1 0.96819 1.000 0.191 -0.942 -0.762 -0.209
2 0.91195 0.191 1.000 0.000 -0.146 -0.912
3 0.96348 -0.942 0.000 1.000 0.718 -0.000
4 0.76233 -0.762 -0.146 0.718 1.000 0.160
5 0.92088 -0.209 -0.912 -0.000 0.160 1.000
**********
** 7 **SET ERR 0.5
**********
**********
** 8 **SET PRINT 1
**********
**********
** 9 **HESSE 2500
**********
COVARIANCE MATRIX CALCULATED SUCCESSFULLY
FCN=16.2872 FROM HESSE STATUS=OK 31 CALLS 222 TOTAL
EDM=3.6099e-06 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER INTERNAL INTERNAL
NO. NAME VALUE ERROR STEP SIZE VALUE
1 b 8.35905e+01 1.38734e+01 1.46135e-05 -4.58641e-01
2 bbar 9.93029e+02 3.14977e+01 1.03809e-05 -3.49712e-01
3 rho 1.27859e+00 1.98839e-01 3.16616e-05 2.82321e-01
4 s 5.54105e+01 1.78190e+01 1.34101e-04 1.08422e-01
5 tau 4.94037e+00 1.71676e-01 1.89816e-05 -2.98215e-02
ERR DEF= 0.5
EXTERNAL ERROR MATRIX. NDIM= 25 NPAR= 5 ERR DEF=0.5
1.932e+02 8.358e+01 -2.623e+00 -1.931e+02 -5.000e-01
8.358e+01 9.928e+02 -2.730e-04 -8.358e+01 -4.939e+00
-2.623e+00 -2.730e-04 4.012e-02 2.623e+00 1.633e-06
-1.931e+02 -8.358e+01 2.623e+00 3.321e+02 5.000e-01
-5.000e-01 -4.939e+00 1.633e-06 5.000e-01 2.955e-02
PARAMETER CORRELATION COEFFICIENTS
NO. GLOBAL 1 2 3 4 5
1 0.96822 1.000 0.191 -0.942 -0.763 -0.209
2 0.91193 0.191 1.000 -0.000 -0.146 -0.912
3 0.96351 -0.942 -0.000 1.000 0.718 0.000
4 0.76256 -0.763 -0.146 0.718 1.000 0.160
5 0.92086 -0.209 -0.912 0.000 0.160 1.000
[#1] INFO:Minization -- RooMinimizer::optimizeConst: deactivating const optimization
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_model_modelData_Profile[s]) Creating instance of MINUIT
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_model_modelData_Profile[s]) determining minimum likelihood for current configurations w.r.t all observable
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_model_modelData_Profile[s]) minimum found at (s=55.4077)
.
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_model_modelData_Profile[s]) Creating instance of MINUIT
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_model_modelData_Profile[s]) determining minimum likelihood for current configurations w.r.t all observable
[#0] ERROR:InputArguments -- RooArgSet::checkForDup: ERROR argument with name s is already in this set
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_model_modelData_Profile[s]) minimum found at (s=55.4105)
..........................................................................................................................................................................................................Profile Likelihood interval on s = [12.1902, 88.6871]
Real time 0:00:00, CP time 0.550
void FourBinInstructional(bool doBayesian = false, bool doFeldmanCousins = false, bool doMCMC = false)
{
wspace->
factory(
"Poisson::on(non[0,1000], sum::splusb(s[40,0,100],b[100,0,300]))");
wspace->
factory(
"Poisson::off(noff[0,5000], prod::taub(b,tau[5,3,7],rho[1,0,2]))");
wspace->
factory(
"Poisson::onbar(nonbar[0,10000], bbar[1000,500,2000])");
wspace->
factory(
"Poisson::offbar(noffbar[0,1000000], prod::lambdaoffbar(bbar, tau))");
wspace->
factory(
"Gaussian::mcCons(rhonom[1.,0,2], rho, sigma[.2])");
wspace->
factory(
"PROD::model(on,off,onbar,offbar,mcCons)");
wspace->
defineSet(
"obs",
"non,noff,nonbar,noffbar,rhonom");
wspace->
factory(
"Uniform::prior_poi({s})");
wspace->
factory(
"Uniform::prior_nuis({b,bbar,tau, rho})");
wspace->
factory(
"PROD::prior(prior_poi,prior_nuis)");
plc.SetConfidenceLevel(0.95);
fc.SetConfidenceLevel(0.95);
fc.FluctuateNumDataEntries(
false);
fc.UseAdaptiveSampling(
true);
if (doFeldmanCousins) {
}
bc.SetConfidenceLevel(0.95);
if (doBayesian && wspace->
set(
"poi")->
getSize() == 1) {
bInt = bc.GetInterval();
} else {
cout << "Bayesian Calc. only supports on parameter of interest" << endl;
}
mc.SetConfidenceLevel(0.95);
mc.SetProposalFunction(*pf);
mc.SetNumBurnInSteps(500);
mc.SetNumIters(50000);
mc.SetLeftSideTailFraction(0.5);
if (doMCMC)
mcInt = mc.GetInterval();
if (doBayesian && doMCMC) {
} else if (doBayesian || doMCMC) {
}
if (doBayesian && wspace->
set(
"poi")->
getSize() == 1) {
bc.SetScanOfPosterior(20);
RooPlot *bplot = bc.GetPosteriorPlot();
}
if (doMCMC) {
if (doBayesian && wspace->
set(
"poi")->
getSize() == 1)
else
mcPlot.Draw();
}
cout <<
"Profile Likelihood interval on s = [" << plInt->
LowerLimit(*wspace->
var(
"s")) <<
", "
if (doBayesian && wspace->
set(
"poi")->
getSize() == 1) {
cout <<
"Bayesian interval on s = [" << bInt->
LowerLimit() <<
", " << bInt->
UpperLimit() <<
"]" << endl;
}
if (doFeldmanCousins) {
cout <<
"Feldman Cousins interval on s = [" << fcInt->
LowerLimit(*wspace->
var(
"s")) <<
", "
}
if (doMCMC) {
cout <<
"MCMC interval on s = [" << mcInt->
LowerLimit(*wspace->
var(
"s")) <<
", "
}
}
static struct mg_connection * fc(struct mg_context *ctx)
RooDataSet * generate(const RooArgSet &whatVars, Int_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none())
See RooAbsPdf::generate(const RooArgSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,...
virtual RooFitResult * fitTo(RooAbsData &data, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Fit PDF to given dataset.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
RooDataSet is a container class to hold unbinned data.
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
static RooMsgService & instance()
Return reference to singleton instance.
void setGlobalKillBelow(RooFit::MsgLevel level)
RooFit::MsgLevel globalKillBelow() const
A RooPlot is a plot frame and a container for graphics objects within that frame.
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
BayesianCalculator is a concrete implementation of IntervalCalculator, providing the computation of a...
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
This class provides simple and straightforward utilities to plot a LikelihoodInterval object.
void Draw(const Option_t *options=0)
draw the likelihood interval or contour for the 1D case a RooPlot is drawn by default of the profiled...
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
Double_t LowerLimit(const RooRealVar ¶m)
return the lower bound of the interval on a given parameter
Double_t UpperLimit(const RooRealVar ¶m)
return the upper bound of the interval on a given parameter
Bayesian Calculator estimating an interval or a credible region using the Markov-Chain Monte Carlo me...
This class provides simple and straightforward utilities to plot a MCMCInterval object.
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
virtual Double_t UpperLimit(RooRealVar ¶m)
get the highest value of param that is within the confidence interval
virtual Double_t LowerLimit(RooRealVar ¶m)
get the lowest value of param that is within the confidence interval
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
virtual void SetPriorPdf(const RooAbsPdf &pdf)
Set the Prior Pdf, add to the the workspace if not already there.
virtual void SetWorkspace(RooWorkspace &ws)
virtual void SetParametersOfInterest(const RooArgSet &set)
Specify parameters of interest.
virtual void SetNuisanceParameters(const RooArgSet &set)
Specify the nuisance parameters (parameters that are not POI).
virtual void SetPdf(const RooAbsPdf &pdf)
Set the Pdf, add to the the workspace if not already there.
PointSetInterval is a concrete implementation of the ConfInterval interface.
Double_t UpperLimit(RooRealVar ¶m)
return upper limit on a given parameter
Double_t LowerLimit(RooRealVar ¶m)
return lower limit on a given parameter
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
ProposalFunction is an interface for all proposal functions that would be used with a Markov Chain Mo...
virtual void SetCovMatrix(const TMatrixDSym &covMatrix)
virtual void SetUpdateProposalParameters(Bool_t updateParams)
virtual ProposalFunction * GetProposalFunction()
virtual void SetVariables(RooArgList &vars)
virtual void SetCacheSize(Int_t size)
SimpleInterval is a concrete implementation of the ConfInterval interface.
virtual Double_t UpperLimit()
virtual Double_t LowerLimit()
The RooWorkspace is a persistable container for RooFit projects.
Bool_t writeToFile(const char *fileName, Bool_t recreate=kTRUE)
Save this current workspace into given file.
Bool_t defineSet(const char *name, const RooArgSet &aset, Bool_t importMissing=kFALSE)
Define a named RooArgSet with given constituents.
RooRealVar * var(const char *name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found.
Bool_t import(const RooAbsArg &arg, const RooCmdArg &arg1=RooCmdArg(), const RooCmdArg &arg2=RooCmdArg(), const RooCmdArg &arg3=RooCmdArg(), const RooCmdArg &arg4=RooCmdArg(), const RooCmdArg &arg5=RooCmdArg(), const RooCmdArg &arg6=RooCmdArg(), const RooCmdArg &arg7=RooCmdArg(), const RooCmdArg &arg8=RooCmdArg(), const RooCmdArg &arg9=RooCmdArg())
Import a RooAbsArg object, e.g.
RooFactoryWSTool & factory()
Return instance to factory tool.
const RooArgSet * set(const char *name)
Return pointer to previously defined named set with given nmame If no such set is found a null pointe...
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
RooCmdArg Save(Bool_t flag=kTRUE)
fit(model, train_loader, val_loader, num_epochs, batch_size, optimizer, criterion, save_best, scheduler)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
Namespace for the RooStats classes.