A hypothesis testing example based on number counting with background uncertainty.
NOTE: This example must be run with the ACLIC (the + option ) due to the new class that is defined.
The basic setup here is that a main measurement has observed x events with an expectation of s+b. One can choose an ad hoc prior for the uncertainty on b, or try to base it on an auxiliary measurement. In this case, the auxiliary measurement (aka control measurement, sideband) is another counting experiment with measurement y and expectation tau*b. With an 'original prior' on b, called \(\eta(b)\) then one can obtain a posterior from the auxiliary measurement \(\pi(b) = \eta(b) * Pois(y|tau*b).\) This is a principled choice for a prior on b in the main measurement of x, which can then be treated in a hybrid Bayesian/Frequentist way. Additionally, one can try to treat the two measurements simultaneously, which is detailed in Part 6 of the tutorial.
This tutorial is related to the FourBin.C tutorial in the modeling, but focuses on hypothesis testing instead of interval estimation.
More background on this 'prototype problem' can be found in the following papers:
[#1] INFO:NumericIntegration -- RooRealIntegral::init([py_X_prior_b_X_px]_Norm[b]_denominator_Int[b]) using numeric integrator RooRombergIntegrator to calculate Int(b)
-----------------------------------------
Part 2
Hybrid p-value from direct integration = 0.000774094
Z_Gamma Significance = 3.1655
-----------------------------------------
Part 3
Z_Bi p-value (analytic): 0.00094165
Z_Bi significance (analytic): 3.10804
Real time 0:00:05, CP time 5.250
[#0] WARNING:InputArguments -- The parameter 'b' with range [0, 300] of the RooLognormal 'lognorm_prior' exceeds the safe range of (0, inf). Advise to limit its range.
[#0] WARNING:InputArguments -- The parameter 'y' with range [0, 500] of the RooLognormal 'lognorm_prior' exceeds the safe range of (0, inf). Advise to limit its range.
[#0] WARNING:Eval -- RooStatsUtils::MakeNuisancePdf - no constraints found on nuisance parameters in the input model
[#0] WARNING:Eval -- RooStatsUtils::MakeNuisancePdf - no constraints found on nuisance parameters in the input model
-----------------------------------------
Part 4
Results HypoTestCalculator_result:
- Null p-value = 0.00115 +/- 0.000239654
- Significance = 3.04848 +/- 0.0626148 sigma
- Number of Alt toys: 1000
- Number of Null toys: 20000
- Test statistic evaluated on data: 150
- CL_b: 0.00115 +/- 0.000239654
- CL_s+b: 0.524 +/- 0.0157932
- CL_s: 455.652 +/- 95.9434
Real time 0:00:02, CP time 2.520
[#0] WARNING:Eval -- RooStatsUtils::MakeNuisancePdf - no constraints found on nuisance parameters in the input model
[#0] WARNING:Eval -- RooStatsUtils::MakeNuisancePdf - no constraints found on nuisance parameters in the input model
-----------------------------------------
Part 5
Results HypoTestCalculator_result:
- Null p-value = 0.0009 +/- 0.000212037
- Significance = 3.12139 +/- 0.0693716 sigma
- Number of Alt toys: 1000
- Number of Null toys: 20000
- Test statistic evaluated on data: 10.8198
- CL_b: 0.0009 +/- 0.000212037
- CL_s+b: 0.535 +/- 0.0157726
- CL_s: 594.444 +/- 141.141
Real time 0:00:02, CP time 2.770
[#0] WARNING:Eval -- RooStatsUtils::MakeNuisancePdf - no constraints found on nuisance parameters in the input model
[#0] WARNING:Eval -- RooStatsUtils::MakeNuisancePdf - no constraints found on nuisance parameters in the input model
-----------------------------------------
Part 6
Results HypoTestCalculator_result:
- Null p-value = 0.000666667 +/- 0.000149021
- Significance = 3.20871 +/- 0.0642752 sigma
- Number of Alt toys: 1000
- Number of Null toys: 30000
- Test statistic evaluated on data: 5.03388
- CL_b: 0.000666667 +/- 0.000149021
- CL_s+b: 0.489 +/- 0.0158076
- CL_s: 733.5 +/- 165.667
Real time 0:00:27, CP time 27.590
public:
BinCountTestStat(void) : fColumnName("tmp") {}
BinCountTestStat(string columnName) : fColumnName(columnName) {}
{
for (
int i = 0; i <
data.numEntries(); i++) {
value +=
data.get(i)->getRealValue(fColumnName.c_str());
}
}
virtual const TString GetVarName()
const {
return fColumnName; }
private:
string fColumnName;
protected:
};
void HybridInstructional()
{
w->factory(
"Poisson::px(x[150,0,500],sum::splusb(s[0,0,100],b[100,0,300]))");
w->factory(
"Poisson::py(y[100,0,500],prod::taub(tau[1.],b))");
w->factory(
"PROD::model(px,py)");
w->factory(
"Uniform::prior_b(b)");
w->factory(
"PROJ::averagedModel(PROD::foo(px|b,py,prior_b),b)");
w->var(
"s")->setVal(50.);
w->var(
"y")->setVal(100);
w->var(
"x")->setVal(150);
std::unique_ptr<RooAbsReal> cdf{
w->pdf(
"averagedModel")->createCdf(*
w->var(
"x"))};
cdf->getVal();
cout << "-----------------------------------------" << endl;
cout << "Part 2" << endl;
cout << "Hybrid p-value from direct integration = " << 1 - cdf->getVal() << endl;
double p_Bi = NumberCountingUtils::BinomialWithTauObsP(150, 100, 1);
double Z_Bi = NumberCountingUtils::BinomialWithTauObsZ(150, 100, 1);
cout << "-----------------------------------------" << endl;
cout << "Part 3" << endl;
std::cout << "Z_Bi p-value (analytic): " << p_Bi << std::endl;
std::cout << "Z_Bi significance (analytic): " << Z_Bi << std::endl;
w->defineSet(
"obs",
"x");
w->defineSet(
"poi",
"s");
data->add(*
w->set(
"obs"));
b_model.SetPdf(*
w->pdf(
"px"));
b_model.SetObservables(*
w->set(
"obs"));
b_model.SetParametersOfInterest(*
w->set(
"poi"));
w->var(
"s")->setVal(0.0);
b_model.SetSnapshot(*
w->set(
"poi"));
sb_model.SetPdf(*
w->pdf(
"px"));
sb_model.SetObservables(*
w->set(
"obs"));
sb_model.SetParametersOfInterest(*
w->set(
"poi"));
w->var(
"s")->setVal(50.0);
sb_model.SetSnapshot(*
w->set(
"poi"));
BinCountTestStat binCount("x");
w->factory(
"Gaussian::gauss_prior(b,y, expr::sqrty('sqrt(y)',y))");
w->factory(
"Lognormal::lognorm_prior(b,y, expr::kappa('1+1./sqrt(y)',y))");
hc1.SetToys(20000, 1000);
hc1.ForcePriorNuisanceAlt(*
w->pdf(
"py"));
hc1.ForcePriorNuisanceNull(*
w->pdf(
"py"));
cout << "-----------------------------------------" << endl;
cout << "Part 4" << endl;
slrts.SetNullParameters(*b_model.GetSnapshot());
slrts.SetAltParameters(*sb_model.GetSnapshot());
hc2.SetToys(20000, 1000);
hc2.ForcePriorNuisanceAlt(*
w->pdf(
"py"));
hc2.ForcePriorNuisanceNull(*
w->pdf(
"py"));
if (pc)
cout << "-----------------------------------------" << endl;
cout << "Part 5" << endl;
w->defineSet(
"obsXY",
"x,y");
w->var(
"x")->setVal(150.);
w->var(
"y")->setVal(100.);
dataXY->
add(*
w->set(
"obsXY"));
b_modelXY.SetPdf(*
w->pdf(
"model"));
b_modelXY.SetObservables(*
w->set(
"obsXY"));
b_modelXY.SetParametersOfInterest(*
w->set(
"poi"));
w->var(
"s")->setVal(0.0);
b_modelXY.SetSnapshot(*
w->set(
"poi"));
sb_modelXY.SetPdf(*
w->pdf(
"model"));
sb_modelXY.SetObservables(*
w->set(
"obsXY"));
sb_modelXY.SetParametersOfInterest(*
w->set(
"poi"));
w->var(
"s")->setVal(50.0);
sb_modelXY.SetSnapshot(*
w->set(
"poi"));
ropl.SetSubtractMLE(false);
w->factory(
"Gamma::gamma_y0(b,sum::temp0(y0,1),1,0)");
w->factory(
"Gaussian::gauss_prior_y0(b,y0, expr::sqrty0('sqrt(y0)',y0))");
hc3.SetToys(30000, 1000);
hc3.ForcePriorNuisanceAlt(*
w->pdf(
"gamma_y0"));
hc3.ForcePriorNuisanceNull(*
w->pdf(
"gamma_y0"));
if (pc)
cout << "-----------------------------------------" << endl;
cout << "Part 6" << endl;
}
#define ClassDef(name, id)
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
Abstract base class for binned and unbinned datasets.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
RooDataSet is a container class to hold unbinned data.
void add(const RooArgSet &row, double weight, double weightError)
Add one ore more rows of data.
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.
static RooPlot * frame(const RooAbsRealLValue &var, double xmin, double xmax, Int_t nBins)
Create a new frame for a given variable in x.
void Draw(Option_t *options=nullptr) override
Draw this plot and all of the elements it contains.
Same purpose as HybridCalculatorOriginal, but different implementation.
This class provides the plots for the result of a study performed with any of the HypoTestCalculatorG...
HypoTestResult is a base class for results from hypothesis tests.
void Print(const Option_t *="") const override
Print out some information about the results Note: use Alt/Null labels for the hypotheses here as the...
MaxLikelihoodEstimateTestStat: TestStatistic that returns maximum likelihood estimate of a specified ...
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
ProfileLikelihoodTestStat is an implementation of the TestStatistic interface that calculates the pro...
Holds configuration options for proof and proof-lite.
TestStatistic that returns the ratio of profiled likelihoods.
void Draw(Option_t *options=nullptr) override
Draw this plot and all of the elements it contains.
TestStatistic class that returns -log(L[null] / L[alt]) where L is the likelihood.
TestStatistic is an interface class to provide a facility for construction test statistics distributi...
ToyMCSampler is an implementation of the TestStatSampler interface.
void SetProofConfig(ProofConfig *pc=nullptr)
calling with argument or nullptr deactivates proof
virtual void SetTestStatistic(TestStatistic *testStatistic, unsigned int i)
Set the TestStatistic (want the argument to be a function of the data & parameter points.
virtual void SetNEventsPerToy(const Int_t nevents)
Forces the generation of exactly n events even for extended PDFs.
Persistable container for RooFit projects.
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
Automatic pad generation by division.
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
void Stop()
Stop the stopwatch.
void Print(Option_t *option="") const override
Print the real and cpu time passed between the start and stop events.
RooCmdArg LineColor(Color_t color)
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.
double PValueToSignificance(double pvalue)
returns one-sided significance corresponding to a p-value