Example demonstrating usage of HybridCalcultor
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.
This example:
- demonstrates the usage of the HybridCalcultor (Part 4-6)
- demonstrates the numerical integration of RooFit (Part 2)
- validates the RooStats against an example with a known analytic answer
- demonstrates usage of different test statistics
- explains subtle choices in the prior used for hybrid methods
- demonstrates usage of different priors for the nuisance parameters
- demonstrates usage of PROOF
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:
Processing /mnt/build/workspace/root-makedoc-v608/rootspi/rdoc/src/v6-08-00-patches/tutorials/roostats/HybridInstructional.C...
public:
BinCountTestStat(void) : fColumnName("tmp") {}
BinCountTestStat(string columnName) : fColumnName(columnName) {}
}
return value;
}
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(
"PROJ::averagedModel(PROD::foo(px|b,py,prior_b),b)") ;
cout << "-----------------------------------------"<<endl;
cout << "Part 2" << endl;
cout <<
"Hybrid p-value from direct integration = " << 1-cdf->
getVal() << endl;
cout << "Z_Gamma Significance = " <<
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;
b_model.SetPdf(*w->
pdf(
"px"));
b_model.SetObservables(*w->
set(
"obs"));
b_model.SetParametersOfInterest(*w->
set(
"poi"));
b_model.SetSnapshot(*w->
set(
"poi"));
sb_model.SetPdf(*w->
pdf(
"px"));
sb_model.SetObservables(*w->
set(
"obs"));
sb_model.SetParametersOfInterest(*w->
set(
"poi"));
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.SetAltParameters(*sb_model.GetSnapshot());
hc2.SetToys(20000,1000);
hc2.ForcePriorNuisanceAlt(*w->
pdf(
"py"));
hc2.ForcePriorNuisanceNull(*w->
pdf(
"py"));
cout << "-----------------------------------------"<<endl;
cout << "Part 5" << endl;
dataXY->
add(*w->
set(
"obsXY"));
b_modelXY.SetPdf(*w->
pdf(
"model"));
b_modelXY.SetObservables(*w->
set(
"obsXY"));
b_modelXY.SetParametersOfInterest(*w->
set(
"poi"));
b_modelXY.SetSnapshot(*w->
set(
"poi"));
sb_modelXY.SetPdf(*w->
pdf(
"model"));
sb_modelXY.SetObservables(*w->
set(
"obsXY"));
sb_modelXY.SetParametersOfInterest(*w->
set(
"poi"));
sb_modelXY.SetSnapshot(*w->
set(
"poi"));
ropl(*b_modelXY.GetPdf(), *sb_modelXY.GetPdf(), sb_modelXY.GetSnapshot());
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"));
cout << "-----------------------------------------"<<endl;
cout << "Part 6" << endl;
}
- Authors
- Kyle Cranmer, Wouter Verkerke, and Sven Kreiss
Definition in file HybridInstructional.C.