diff --git a/roofit/roofitcore/test/CMakeLists.txt b/roofit/roofitcore/test/CMakeLists.txt index 810eb310a7fb4..34d6338740196 100644 --- a/roofit/roofitcore/test/CMakeLists.txt +++ b/roofit/roofitcore/test/CMakeLists.txt @@ -94,8 +94,7 @@ ROOT_ADD_GTEST(testRooTruthModel testRooTruthModel.cxx LIBRARIES RooFitCore RooF COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/rooAbsAnaConvPdf_classV3.root) if (roofit_multiprocess) - ROOT_ADD_GTEST(testTestStatisticsPlot TestStatistics/testPlot.cxx LIBRARIES RooFitMultiProcess RooFitCore RooFit - COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/TestStatistics/TestStatistics_ref.root) + ROOT_ADD_GTEST(testTestStatisticsPlot TestStatistics/testPlot.cxx LIBRARIES RooFitMultiProcess RooFitCore RooFit) ROOT_ADD_GTEST(testLikelihoodGradientJob TestStatistics/testLikelihoodGradientJob.cxx LIBRARIES RooFitMultiProcess RooFitCore m ROOT::TestSupport) target_include_directories(testLikelihoodGradientJob PRIVATE ${RooFitCore_MultiProcess_TestStatistics_INCLUDE_DIR}) ROOT_ADD_GTEST(testLikelihoodJob TestStatistics/testLikelihoodJob.cxx LIBRARIES RooFitMultiProcess RooFitCore m) diff --git a/roofit/roofitcore/test/TestStatistics/TestStatistics_ref.root b/roofit/roofitcore/test/TestStatistics/TestStatistics_ref.root deleted file mode 100644 index 6b8c89029b9ea..0000000000000 Binary files a/roofit/roofitcore/test/TestStatistics/TestStatistics_ref.root and /dev/null differ diff --git a/roofit/roofitcore/test/TestStatistics/testPlot.cxx b/roofit/roofitcore/test/TestStatistics/testPlot.cxx index fbedc0968468b..8a3ee13948659 100644 --- a/roofit/roofitcore/test/TestStatistics/testPlot.cxx +++ b/roofit/roofitcore/test/TestStatistics/testPlot.cxx @@ -11,83 +11,99 @@ */ #include +#include #include -#include -#include -#include +#include #include #include #include #include -#include #include -#include - #include +#include #include using namespace RooFit; -class TestRooRealLPlot : public RooUnitTest { -public: - TestRooRealLPlot(TFile &refFile, bool writeRef, int verbose) - : RooUnitTest("Plotting and minimization with RooFit::TestStatistics", &refFile, writeRef, verbose){}; - bool testCode() override - { - - // C r e a t e m o d e l a n d d a t a - // --------------------------------------- - // Constructing a workspace with pdf and dataset - RooWorkspace w("w"); - w.factory("expr::Nexp('mu*S+B',mu[1,-1,10],S[10],B[20])"); - w.factory("Poisson::model(Nobs[0,100],Nexp)"); - w.var("Nobs")->setBins(4); - RooDataSet d("d", "d", *w.var("Nobs")); - w.var("Nobs")->setVal(25); - d.add(*w.var("Nobs")); - - // P e r f o r m a p a r a l l e l l i k e l i h o o d m i n i m i z a t i o n - // -------------------------------------------------------------------------------- - - // Creating a RooAbsL likelihood - std::unique_ptr likelihood{w.pdf("model")->createNLL(d, ModularL(true))}; - - // Creating a minimizer and explicitly setting type of parallelization - std::size_t nWorkers = 1; - RooMinimizer::Config cfg; - cfg.parallelize = nWorkers; - cfg.enableParallelDescent = false; - cfg.enableParallelGradient = true; - RooMinimizer m(*likelihood, cfg); - - // Minimize - m.migrad(); - - // C o n v e r t t o R o o R e a l L a n d p l o t - // --------------------------------------------------- - RooPlot *xframe = w.var("mu")->frame(-1, 10); - likelihood->plotOn(xframe, RooFit::Precision(1)); - - // --- Post processing for RooUnitTest --- - regPlot(xframe, "TestRooRealLPlot_plot"); - - return true; - } -}; - +/// Plotting and minimization with RooFit::TestStatistics: minimize a +/// RooFit::TestStatistics likelihood with the parallel gradient and plot it as +/// a function of the parameter. The results are validated against analytic +/// expectations and against direct evaluations of the likelihood, instead of +/// the RooUnitTest reference file that was used before. TEST(TestStatisticsPlot, RooRealL) { - // Run the RooUnitTest and assert that it succeeds with gtest - - RooUnitTest::setMemDir(gDirectory); - - gErrorIgnoreLevel = kWarning; - - TFile fref("TestStatistics_ref.root"); + RooHelpers::LocalChangeMsgLevel changeMsgLvl{RooFit::WARNING}; + + // C r e a t e m o d e l a n d d a t a + // --------------------------------------- + // Constructing a workspace with pdf and dataset + RooWorkspace w("w"); + w.factory("expr::Nexp('mu*S+B',mu[1,-1,10],S[10],B[20])"); + w.factory("Poisson::model(Nobs[0,100],Nexp)"); + w.var("Nobs")->setBins(4); + RooDataSet d("d", "d", *w.var("Nobs")); + w.var("Nobs")->setVal(25); + d.add(*w.var("Nobs")); + + RooRealVar &mu = *w.var("mu"); + + // P e r f o r m a p a r a l l e l l i k e l i h o o d m i n i m i z a t i o n + // -------------------------------------------------------------------------------- + + // Creating a RooAbsL likelihood + std::unique_ptr likelihood{w.pdf("model")->createNLL(d, ModularL(true))}; + + // Creating a minimizer and explicitly setting type of parallelization + std::size_t nWorkers = 1; + RooMinimizer::Config cfg; + cfg.parallelize = nWorkers; + cfg.enableParallelDescent = false; + cfg.enableParallelGradient = true; + RooMinimizer m(*likelihood, cfg); + + // Minimize + m.setPrintLevel(-1); + m.migrad(); + + // The analytic maximum likelihood estimate is at Nexp = Nobs, i.e. + // mu = (Nobs - B) / S. The tolerance is at the scale of the Minuit + // convergence criterion, given that sigma(mu) = sqrt(Nobs) / S = 0.5. + EXPECT_NEAR(mu.getVal(), 0.5, 0.05 * 0.5); + + // C o n v e r t t o R o o R e a l L a n d p l o t + // --------------------------------------------------- + std::unique_ptr xframe{mu.frame(-1, 10)}; + likelihood->plotOn(xframe.get(), RooFit::Precision(1)); + RooCurve *curve = xframe->getCurve(); + ASSERT_NE(curve, nullptr); + ASSERT_GT(curve->GetN(), 1); + + // Every point of the plotted curve must match a direct evaluation of the + // likelihood at that parameter value + for (int i = 0; i < curve->GetN(); ++i) { + const double muVal = curve->GetPointX(i); + if (muVal < mu.getMin() || muVal > mu.getMax()) + continue; + mu.setVal(muVal); + const double directVal = likelihood->getVal(); + EXPECT_NEAR(curve->GetPointY(i), directVal, 1e-6 * std::max(1.0, std::abs(directVal))) + << "curve point " << i << " at mu = " << muVal; + } - TestRooRealLPlot plotTest{fref, false, 0}; - bool result = plotTest.runTest(); - ASSERT_TRUE(result); + // The likelihood must also match the analytically known Poisson -log L, + // which is defined up to a mu-independent constant. The comparison is + // restricted to moderate Nexp values, where the truncation of the Poisson + // normalization to the observable range is negligible. + auto analyticNll = [](double muVal) { + const double nexp = 10 * muVal + 20; + return nexp - 25 * std::log(nexp); + }; + mu.setVal(0.5); + const double nllOffset = likelihood->getVal() - analyticNll(0.5); + for (double muVal : {0.0, 1.0, 2.0, 3.0}) { + mu.setVal(muVal); + EXPECT_NEAR(likelihood->getVal(), analyticNll(muVal) + nllOffset, 1e-4) << "at mu = " << muVal; + } } diff --git a/roofit/roostats/test/CMakeLists.txt b/roofit/roostats/test/CMakeLists.txt index e14c02d8a35e6..1896a906e15f9 100644 --- a/roofit/roostats/test/CMakeLists.txt +++ b/roofit/roostats/test/CMakeLists.txt @@ -6,20 +6,19 @@ ROOT_ADD_GTEST(testHypoTestInvResult testHypoTestInvResult.cxx ROOT_ADD_GTEST(testSPlot testSPlot.cxx LIBRARIES RooStats) #--stressRooStats---------------------------------------------------------------------------------- -ROOT_EXECUTABLE(stressRooStats stressRooStats.cxx LIBRARIES RooStats Gpad Net) +# Googletest version of the old RooStats S.T.R.E.S.S. suite, parameterized +# over the RooFit evaluation backends. Like in the original suite, the tests +# are also run with the Minuit2 minimizer (the default minimizer of the suite +# is Minuit2). +ROOT_ADD_GTEST(stressRooStats stressRooStats.cxx LIBRARIES RooStats LABELS longtest TIMEOUT 3600) if(mathmore) target_compile_definitions(stressRooStats PRIVATE ROOFITMORE) endif() - -configure_file(stressRooStats_ref.root stressRooStats_ref.root COPYONLY) -if(roofit_legacy_eval_backend) - ROOT_ADD_TEST(test-stressroostats-legacy COMMAND stressRooStats -b legacy FAILREGEX "FAILED|Error in" LABELS longtest) -endif() -ROOT_ADD_TEST(test-stressroostats-cpu COMMAND stressRooStats -b cpu FAILREGEX "FAILED|Error in" LABELS longtest) if(cuda) - ROOT_ADD_TEST(test-stressroostats-cuda COMMAND stressRooStats -b cuda FAILREGEX "FAILED|Error in" LABELS longtest RESOURCE_LOCK GPU) -endif() -if(roofit_legacy_eval_backend) - ROOT_ADD_TEST(test-stressroostats-legacy-minuit2 COMMAND stressRooStats -minim Minuit2 -b legacy FAILREGEX "FAILED|Error in" LABELS longtest) + set_tests_properties(gtest-roofit-roostats-stressRooStats PROPERTIES RESOURCE_LOCK GPU) endif() -ROOT_ADD_TEST(test-stressroostats-cpu-minuit2 COMMAND stressRooStats -minim Minuit2 -b cpu FAILREGEX "FAILED|Error in" LABELS longtest) +# Like the non-default-minimizer variant in the original suite, the Minuit +# variant only runs the legacy and cpu backends (that is also why it doesn't +# need to lock the GPU resource). +ROOT_ADD_TEST(test-stressroostats-minuit COMMAND stressRooStats --gtest_filter=-*EvalBackendcuda* + ENVIRONMENT STRESSROOSTATS_MINIMIZER=Minuit LABELS longtest TIMEOUT 3600) diff --git a/roofit/roostats/test/stressRooStats.cxx b/roofit/roostats/test/stressRooStats.cxx index 916fb9348eea7..f72b1bf2889a3 100644 --- a/roofit/roostats/test/stressRooStats.cxx +++ b/roofit/roostats/test/stressRooStats.cxx @@ -1,421 +1,1267 @@ -// @(#)root/roofitcore:$name: $:$id$ -// Authors: Wouter Verkerke November 2007 +// Authors: Ioan Gabriel Bucur, Lorenzo Moneta, Wouter Verkerke +// +// Googletest version of the old RooStats S.T.R.E.S.S. suite. Each of the 48 +// tests of the original suite is translated to a (parameterized) gtest test +// case with the same models, calculator configurations, random seeds and +// tolerances. +// +// The original suite compared results with references stored in +// stressRooStats_ref.root. Depending on the test, these references were either +// computed analytically, taken from a publication, or produced by the tested +// calculator itself at reference-writing time (pure regression tests). In this +// version, analytic references are computed inline, and published or +// regression reference values are hardcoded (the latter extracted from the +// last stressRooStats_ref.root). +// +// The tests are parameterized over the RooFit evaluation backends, matching +// the backend coverage of the removed stressRooStats invocations. Like the +// old "-minim" command line option, the minimizer can be chosen via the +// STRESSROOSTATS_MINIMIZER environment variable (default is Minuit2). + +#include "../../roofitcore/test/gtest_wrapper.h" + +// Global functions that build the more complex RooStats models +#include "stressRooStats_models.h" + +// RooStats headers +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include -// C/C++ headers -#include -#include -#include -#include -#include - -// Math headers -#include "Math/MinimizerOptions.h" +// RooFit headers +#include +#include +#include +#include +#include +#include +#include +#include +#include // ROOT headers -#include "TSystem.h" -#include "TString.h" -#include "TStopwatch.h" -#include "TROOT.h" -#include "TLine.h" -#include "TFile.h" -#include "TClass.h" -#include "TF1.h" -#include "TBenchmark.h" - -// RooFit headers -#include "RooGlobalFunc.h" -#include "RooNumIntConfig.h" -#include "RooMsgService.h" -#include "RooResolutionModel.h" -#include "RooRandom.h" +#include +#include +#include +#include +#include -// Tests file -#include "stressRooStats_tests.h" +#include +#include +#include +#include +#include -using std::string, std::list, std::setw, std::setfill, std::left; +using namespace ROOT::Math; using namespace RooFit; - -//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*// -// // -// RooStats Unit Test S.T.R.E.S.S. Suite // -// Authors: Ioan Gabriel Bucur, Lorenzo Moneta, Wouter Verkerke // -// // -//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*// - -//////////////////////////////////////////////////////////////////////////////// -/// Print test program number and its title - -void StatusPrint(const int id, const TString &title, const int status, const int lineWidth) +using namespace RooStats; + +namespace { + +enum ECalculatorType { kAsymptotic = 0, kFrequentist = 1, kHybrid = 2 }; +enum ETestStatType { + kSimpleLR = 0, + kRatioLR = 1, + kProfileLR = 2, + kProfileLROneSided = 3, + kProfileLROneSidedDiscovery = 4 +}; + +// Confidence levels corresponding to one, two and three Gaussian sigmas +const double kCL1Sigma = 2 * normal_cdf(1) - 1; +const double kCL2Sigma = 2 * normal_cdf(2) - 1; +const double kCL3Sigma = 2 * normal_cdf(3) - 1; + +// Value test tolerance of the original RooUnitTest +const double kVTol = 1e-3; + +std::unique_ptr buildHypoTestCalculator(const ECalculatorType calculatorType, + RooAbsData &data, const ModelConfig &nullModel, + const ModelConfig &altModel, const UInt_t toysNull, + const UInt_t toysAlt) { - TString header = TString::Format("Test %d : %s ", id, title.Data()); - std::cout << left << setw(lineWidth) << setfill('.') << header << " " - << (status > 0 ? "OK" : (status < 0 ? "SKIPPED" : "FAILED")) << std::endl; + if (calculatorType == kAsymptotic) { + return std::make_unique(data, altModel, nullModel); + } else if (calculatorType == kFrequentist) { + auto fc = std::make_unique(data, altModel, nullModel); + // set toys for speedup + fc->SetToys(toysNull, toysAlt); + return fc; + } + // kHybrid + auto hc = std::make_unique(data, altModel, nullModel); + // set toys for speedup + hc->SetToys(toysNull, toysAlt); + return hc; } -//////////////////////////////////////////////////////////////////////////////// -/// width of lines when printing test results - -int stressRooStats(const char *refFile, bool writeRef, int verbose, bool allTests, bool oneTest, int testNumber, - bool dryRun, bool doDump, bool doTreeStore) +std::unique_ptr +buildTestStatistic(const ETestStatType testStatType, const ModelConfig &nullModel, const ModelConfig &altModel) { - const int lineWidth = 120; - - // Save memory directory location - auto memDir = gDirectory; - RooUnitTest::setMemDir(gDirectory); + if (testStatType == kSimpleLR) { + auto slrts = std::make_unique(*nullModel.GetPdf(), *altModel.GetPdf()); + if (nullModel.GetSnapshot()) { + RooArgSet nullParams(*nullModel.GetSnapshot()); + if (nullModel.GetNuisanceParameters()) + nullParams.add(*nullModel.GetNuisanceParameters()); + slrts->SetNullParameters(nullParams); + } + if (altModel.GetSnapshot()) { + RooArgSet altParams(*altModel.GetSnapshot()); + if (altModel.GetNuisanceParameters()) + altParams.add(*altModel.GetNuisanceParameters()); + slrts->SetAltParameters(altParams); + } + slrts->SetAlwaysReuseNLL(true); + return slrts; + } else if (testStatType == kRatioLR) { + auto roplts = std::make_unique(*nullModel.GetPdf(), *altModel.GetPdf(), + altModel.GetSnapshot()); + roplts->SetSubtractMLE(false); + roplts->SetAlwaysReuseNLL(true); + return roplts; + } + // kProfileLR, kProfileLROneSided and kProfileLROneSidedDiscovery + auto plts = std::make_unique(*nullModel.GetPdf()); + if (testStatType == kProfileLROneSided) { + plts->SetOneSided(true); + } else if (testStatType == kProfileLROneSidedDiscovery) { + plts->SetOneSidedDiscovery(true); + } + plts->SetAlwaysReuseNLL(true); + return plts; +} - if (doTreeStore) { - RooAbsData::setDefaultStorageType(RooAbsData::Tree); +/// Create the Poisson product model in the workspace and add the observed +/// values to its data set. Returns the "S+B" model configuration (or nullptr +/// if the workspace content is unexpected, to be caught by the caller). +ModelConfig *setupPoissonProductModel(RooWorkspace &ws, int obsValueX, int obsValueY) +{ + buildPoissonProductModel(&ws); + auto model = dynamic_cast(ws.obj("S+B")); + if (model) { + ws.var("x")->setVal(obsValueX); + ws.var("y")->setVal(obsValueY); + ws.data("data")->add(*model->GetObservables()); } + return model; +} - TFile *fref = nullptr; - if (!dryRun) { - if (TString(refFile).Contains("http:")) { - if (writeRef) { - std::cout << "stressRooStats ERROR: reference file must be local file in writing mode" << std::endl; - return -1; +/// Global setup that mirrors the environment of the original stressRooStats +/// executable: minimizer selection and RooFit message streams silenced below +/// the ERROR level. +class StressRooStatsEnvironment : public ::testing::Environment { +public: + void SetUp() override + { + const char *minimizer = gSystem->Getenv("STRESSROOSTATS_MINIMIZER"); + ROOT::Math::MinimizerOptions::SetDefaultMinimizer(minimizer ? minimizer : "Minuit2"); + + // Disable RooFit messages below the ERROR level, but keep a dedicated + // error stream active so that problems remain visible in the test log. + auto &msgSvc = RooMsgService::instance(); + msgSvc.setSilentMode(true); + for (int i = 0; i < msgSvc.numStreams(); ++i) { + if (msgSvc.getStream(i).minLevel < RooFit::ERROR) { + msgSvc.setStreamStatus(i, false); } - fref = TFile::Open(refFile); - } else { - fref = new TFile(refFile, writeRef ? "RECREATE" : ""); - } - if (fref->IsZombie()) { - std::cout << "stressRooStats ERROR: cannot open reference file " << refFile << std::endl; - return -1; } - } + msgSvc.addStream(RooFit::ERROR); - if (dryRun) { - // Preload singletons here so they don't show up in trace accounting - RooNumIntConfig::defaultConfig(); - RooResolutionModel::identity(); - } + AsymptoticCalculator::SetPrintLevel(0); - // Add dedicated logging stream for errors that will remain active in silent mode - RooMsgService::instance().addStream(RooFit::ERROR); - - std::cout << left << setw(lineWidth) << setfill('*') << "" << std::endl; - std::cout << "*" << setw(lineWidth - 2) << setfill(' ') << " RooStats S.T.R.E.S.S. suite " - << "*" << std::endl; - std::cout << setw(lineWidth) << setfill('*') << "" << std::endl; - std::cout << setw(lineWidth) << setfill('*') << "" << std::endl; - - TStopwatch timer; - timer.Start(); - - list testList; - - // 1-5 TEST PLC CONFINT SIMPLE GAUSSIAN : Confidence Level range is (0,1) - testList.push_back(new TestProfileLikelihoodCalculator1(fref, writeRef, verbose, 0.99999)); // boundary case CL -> 1 - testList.push_back( - new TestProfileLikelihoodCalculator1(fref, writeRef, verbose, 2 * ROOT::Math::normal_cdf(3) - 1)); // 3 sigma - testList.push_back( - new TestProfileLikelihoodCalculator1(fref, writeRef, verbose, 2 * ROOT::Math::normal_cdf(2) - 1)); // 2 sigma - testList.push_back( - new TestProfileLikelihoodCalculator1(fref, writeRef, verbose, 2 * ROOT::Math::normal_cdf(1) - 1)); // 1 sigma - testList.push_back(new TestProfileLikelihoodCalculator1(fref, writeRef, verbose, 0.00001)); // boundary case CL -> 0 - - // 6-10 TEST PLC CONFINT SIMPLE POISSON : Observed value range is [0,1000] - testList.push_back(new TestProfileLikelihoodCalculator2(fref, writeRef, verbose, 0)); // boundary Poisson value (0) - testList.push_back(new TestProfileLikelihoodCalculator2(fref, writeRef, verbose, 1)); - testList.push_back(new TestProfileLikelihoodCalculator2(fref, writeRef, verbose, 5)); - testList.push_back(new TestProfileLikelihoodCalculator2(fref, writeRef, verbose, 100)); - testList.push_back(new TestProfileLikelihoodCalculator2(fref, writeRef, verbose, 800)); // boundary Poisson value - - // 11-13 TEST PLC CONFINT PRODUCT POISSON : Observed value range is [0,30] for x=s+b and [0,80] for y=2*s*1.2^beta - testList.push_back(new TestProfileLikelihoodCalculator3(fref, writeRef, verbose, 10, 30)); - testList.push_back(new TestProfileLikelihoodCalculator3(fref, writeRef, verbose, 20, 25)); - testList.push_back( - new TestProfileLikelihoodCalculator3(fref, writeRef, verbose, 15, 20, 2 * ROOT::Math::normal_cdf(2) - 1)); - - // 14 TEST PLC HYPOTEST ON/OFF MODEL - testList.push_back(new TestProfileLikelihoodCalculator4(fref, writeRef, verbose)); - - // 15-18 TEST BC CONFINT CENTRAL SIMPLE POISSON : Observed value range is [0,100] - testList.push_back(new TestBayesianCalculator1(fref, writeRef, verbose, 1)); - testList.push_back(new TestBayesianCalculator1(fref, writeRef, verbose, 3)); - testList.push_back(new TestBayesianCalculator1(fref, writeRef, verbose, 10)); - testList.push_back(new TestBayesianCalculator1(fref, writeRef, verbose, 50)); - - // 19 TEST BC CONFINT SHORTEST SIMPLE POISSON - testList.push_back(new TestBayesianCalculator2(fref, writeRef, verbose)); - - // 20-22 TEST BC CONFINT CENTRAL PRODUCT POISSON : Observed value range is [0,30] for x=s+b and [0,80] for - // y=2*s*1.2^beta - testList.push_back(new TestBayesianCalculator3(fref, writeRef, verbose, 10, 30)); - testList.push_back(new TestBayesianCalculator3(fref, writeRef, verbose, 20, 25)); - testList.push_back(new TestBayesianCalculator3(fref, writeRef, verbose, 15, 20, 2 * ROOT::Math::normal_cdf(2) - 1)); - - // 23-25 TEST MCMCC CONFINT PRODUCT POISSON : Observed value range is [0,30] for x=s+b and [0,80] for y=2*s*1.2^beta - testList.push_back(new TestMCMCCalculator(fref, writeRef, verbose, 10, 30)); - testList.push_back(new TestMCMCCalculator(fref, writeRef, verbose, 20, 25)); - testList.push_back(new TestMCMCCalculator(fref, writeRef, verbose, 15, 20, 2 * ROOT::Math::normal_cdf(2) - 1)); - - // 26 TEST ZBI SIGNIFICANCE - testList.push_back(new TestZBi(fref, writeRef, verbose)); - - // 27-31 TEST PLC VS AC SIGNIFICANCE : Observed value range is [0,300] for on source and [0,1100] for off-source; tau - // has the range [0.1,5.0] - testList.push_back(new TestHypoTestCalculator1(fref, writeRef, verbose, 150, 100, 1.0)); - testList.push_back(new TestHypoTestCalculator1(fref, writeRef, verbose, 200, 100, 1.0)); - testList.push_back(new TestHypoTestCalculator1(fref, writeRef, verbose, 105, 100, 1.0)); - testList.push_back(new TestHypoTestCalculator1(fref, writeRef, verbose, 150, 10, 0.1)); - testList.push_back(new TestHypoTestCalculator1(fref, writeRef, verbose, 150, 400, 4.0)); - - // 32-36 TEST HTC SIGNIFICANCE - testList.push_back(new TestHypoTestCalculator2(fref, writeRef, verbose, kAsymptotic)); - testList.push_back(new TestHypoTestCalculator2(fref, writeRef, verbose, kFrequentist, kSimpleLR)); - testList.push_back(new TestHypoTestCalculator2(fref, writeRef, verbose, kFrequentist, kRatioLR)); - testList.push_back(new TestHypoTestCalculator2(fref, writeRef, verbose, kFrequentist, kProfileLROneSidedDiscovery)); - testList.push_back(new TestHypoTestCalculator2(fref, writeRef, verbose, kHybrid, kProfileLROneSidedDiscovery)); - - // 37-43 TEST HTI PRODUCT POISSON : Observed value range is [0,30] for x=s+b and [0,80] for y=2*s*1.2^beta - testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, kAsymptotic, kProfileLR, 10, 30)); - testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, kAsymptotic, kProfileLR, 20, 25)); - testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, kAsymptotic, kProfileLR, 15, 20)); - testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, kFrequentist, kProfileLR, 10, 30)); - testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, kFrequentist, kProfileLR, 20, 25)); - testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, kFrequentist, kProfileLR, 15, 20)); - testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, kHybrid, kProfileLR, 10, 30)); - - // 44-48 TEST HTI S+B+E POISSON : Observed value range is [0,50] for x = e*s+b - testList.push_back(new TestHypoTestInverter2(fref, writeRef, verbose, kAsymptotic, kProfileLROneSided, 10, 0.95)); - testList.push_back(new TestHypoTestInverter2(fref, writeRef, verbose, kAsymptotic, kProfileLROneSided, 20)); - // testList.push_back(new TestHypoTestInverter2(fref, writeRef, verbose, kFrequentist, kSimpleLR, 10)); - // testList.push_back(new TestHypoTestInverter2(fref, writeRef, verbose, kFrequentist, kSimpleLR, 20)); - testList.push_back(new TestHypoTestInverter2(fref, writeRef, verbose, kFrequentist, kRatioLR, 10, 0.95)); - testList.push_back(new TestHypoTestInverter2(fref, writeRef, verbose, kFrequentist, kProfileLROneSided, 10, 0.95)); - testList.push_back(new TestHypoTestInverter2(fref, writeRef, verbose, kHybrid, kSimpleLR, 10, 0.95)); - - TString suiteType = TString::Format( - " Starting S.T.R.E.S.S. %s", - allTests ? "full suite" : (oneTest ? TString::Format("test %d", testNumber).Data() : "basic suite")); - - std::cout << "*" << setw(lineWidth - 3) << setfill(' ') << suiteType << " *" << std::endl; - std::cout << setw(lineWidth) << setfill('*') << "" << std::endl; - - if (doDump) { - TFile fdbg("stressRooStats_DEBUG.root", "RECREATE"); + // NOTE: RooIntegrator1D is too slow and gives poor results +#ifdef ROOFITMORE + RooAbsReal::defaultIntegratorConfig()->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D"); +#endif } +}; + +[[maybe_unused]] const auto gStressRooStatsEnv = ::testing::AddGlobalTestEnvironment(new StressRooStatsEnvironment); + +// Reset random generator seeds to make results independent of test ordering, +// like RooUnitTest::runTest() did in the original suite. +void setUpStressRooStatsTest() +{ + gRandom->SetSeed(12345); + RooRandom::randomGenerator()->SetSeed(12345); + RooMsgService::instance().clearErrorCount(); +} - gBenchmark->Start("stressRooStats"); +// The original suite failed a test if RooFit ERROR messages were logged. +void tearDownStressRooStatsTest() +{ + EXPECT_EQ(RooMsgService::instance().errorCount(), 0) << "RooFit ERROR messages were logged during the test"; +} - int nFailed = 0; +/// Fixture for tests that only take the evaluation backend as parameter. The +/// RooStats calculators create their likelihoods internally, so the backend +/// parameter has to be applied via the global default backend (restored again +/// in TearDown to not leak it into other tests). +class StressRooStatsBackendTest : public ::testing::TestWithParam> { +protected: + void SetUp() override { - int i; - list::iterator iter; - - if (oneTest && (testNumber <= 0 || (UInt_t)testNumber > testList.size())) { - std::cout << "Tests are numbered from 1 to " << testList.size() << std::endl; - } else { - for (iter = testList.begin(), i = 1; iter != testList.end(); iter++, i++) { - if (!oneTest || testNumber == i) { - if (doDump) { - (*iter)->setDebug(true); - } - int status = (*iter)->isTestAvailable() ? (*iter)->runTest() : -1; - StatusPrint(i, (*iter)->GetName(), status, lineWidth); - if (!status) - nFailed++; // do not count the skipped tests - } - delete *iter; - } - } + _prevBackend = RooFit::EvalBackend::defaultValue(); + RooFit::EvalBackend::defaultValue() = std::get<0>(GetParam()).value(); + setUpStressRooStatsTest(); + } + void TearDown() override + { + tearDownStressRooStatsTest(); + RooFit::EvalBackend::defaultValue() = _prevBackend; } - gBenchmark->Stop("stressRooStats"); - - // Print table with results - bool UNIX = strcmp(gSystem->GetName(), "Unix") == 0; - std::cout << setw(lineWidth) << setfill('*') << "" << std::endl; - if (UNIX) { - TString sp = gSystem->GetFromPipe("uname -a"); - std::cout << "* SYS: " << sp << std::endl; - if (strstr(gSystem->GetBuildNode(), "Darwin")) { - sp = gSystem->GetFromPipe("sw_vers -productVersion"); - sp += " Mac OS X "; - std::cout << "* SYS: " << sp << std::endl; - } - } else { - const Char_t *os = gSystem->Getenv("OS"); - if (!os) { - std::cout << "* SYS: Windows 95" << std::endl; - } else { - std::cout << "* SYS: " << os << " " << gSystem->Getenv("PROCESSOR_IDENTIFIER") << std::endl; - } +private: + RooFit::EvalBackend::Value _prevBackend = RooFit::EvalBackend::Value::Cpu; +}; + +/// Fixture for tests that are parameterized over the evaluation backend and a +/// test-specific parameter set (see StressRooStatsBackendTest). +template +class StressRooStatsParamTest : public ::testing::TestWithParam> { +protected: + void SetUp() override + { + _prevBackend = RooFit::EvalBackend::defaultValue(); + RooFit::EvalBackend::defaultValue() = std::get<0>(this->GetParam()).value(); + setUpStressRooStatsTest(); + } + void TearDown() override + { + tearDownStressRooStatsTest(); + RooFit::EvalBackend::defaultValue() = _prevBackend; } - std::cout << setw(lineWidth) << setfill('*') << "" << std::endl; - gBenchmark->Print("stressRooStats"); -#ifdef __CLING__ - Double_t reftime = 186.34; // pcbrun4 interpreted -#else - Double_t reftime = 93.59; // pcbrun4 compiled -#endif - const Double_t rootmarks = 860 * reftime / gBenchmark->GetCpuTime("stressRooStats"); + const ParamType ¶m() const { return std::get<1>(this->GetParam()); } + +private: + RooFit::EvalBackend::Value _prevBackend = RooFit::EvalBackend::Value::Cpu; +}; + +std::string backendName(const ::testing::TestParamInfo> &info) +{ + return "EvalBackend" + std::get<0>(info.param).name(); +} - std::cout << setw(lineWidth) << setfill('*') << "" << std::endl; - std::cout << TString::Format("* ROOTMARKS = %6.1f * Root %-8s %d/%d", rootmarks, gROOT->GetVersion(), - gROOT->GetVersionDate(), gROOT->GetVersionTime()) - << std::endl; - std::cout << setw(lineWidth) << setfill('*') << "" << std::endl; +template +std::string backendParamName(const ::testing::TestParamInfo> &info) +{ + return "EvalBackend" + std::get<0>(info.param).name() + "_" + std::get<1>(info.param).name; +} - // NOTE: The function TStopwatch::CpuTime() calls Tstopwatch::Stop(), so you do not need to stop the timer - // separately. - std::cout << "Time at the end of job = " << timer.CpuTime() << " seconds" << std::endl; +} // namespace + +/////////////////////////////////////////////////////////////////////////////// +// +// PART ONE: PROFILE LIKELIHOOD CALCULATOR UNIT TESTS +// +/////////////////////////////////////////////////////////////////////////////// + +/////////////////////////////////////////////////////////////////////////////// +// +// PROFILE LIKELIHOOD CALCULATOR - LIKELIHOOD INTERVAL - GAUSSIAN DISTRIBUTION +// +// Test the likelihood interval computed by the profile likelihood calculator +// on a Gaussian distribution. Reference interval limits are computed via +// analytic methods: solve equation 2*(ln(LL(xMax))-ln(LL(x)) = q, where q = +// normal_quantile_c(testSize/2, 1). In the case of a Gaussian distribution, +// the interval limits are equal to: +// mean +- normal_quantile_c(testSize/2, sigma/sqrt(N)). +// +// ModelConfig (implicit) : +// Observable -> x +// Parameter of Interest -> mean +// Nuisance parameter (constant!) -> sigma +// +// Original tests 1-5 (TestProfileLikelihoodCalculator1), with the confidence +// level probing the boundaries of the (0,1) range. +// +/////////////////////////////////////////////////////////////////////////////// + +struct PlcGaussianParams { + std::string name; + double confidenceLevel; +}; + +class PlcGaussianInterval : public StressRooStatsParamTest {}; + +TEST_P(PlcGaussianInterval, CompareWithAnalyticInterval) +{ + const double confidenceLevel = param().confidenceLevel; + const int N = 10; // number of observations + + // Create Gaussian model and generate a data set + RooWorkspace ws{"w"}; + ws.factory("Gaussian::gauss(x[-5,5], mean[0,-5,5], sigma[1])"); + std::unique_ptr data{ws.pdf("gauss")->generate(*ws.var("x"), N)}; + + // Reference likelihood interval limits computed via analytic methods + const double estMean = data->mean(*ws.var("x")); + const double intervalHalfWidth = + normal_quantile_c((1.0 - confidenceLevel) / 2.0, ws.var("sigma")->getValV() / std::sqrt((double)N)); + + // Calculate likelihood interval using the ProfileLikelihoodCalculator + ProfileLikelihoodCalculator plc{*data, *ws.pdf("gauss"), *ws.var("mean")}; + plc.SetConfidenceLevel(confidenceLevel); + std::unique_ptr interval{plc.GetInterval()}; + + EXPECT_NEAR(interval->LowerLimit(*ws.var("mean")), estMean - intervalHalfWidth, kVTol); + EXPECT_NEAR(interval->UpperLimit(*ws.var("mean")), estMean + intervalHalfWidth, kVTol); +} - if (fref) { - fref->Close(); - delete fref; +INSTANTIATE_TEST_SUITE_P(RooStats, PlcGaussianInterval, + ::testing::Combine(::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values( + PlcGaussianParams{"CLNearOne", 0.99999}, // boundary case CL -> 1 + PlcGaussianParams{"CL3Sigma", kCL3Sigma}, + PlcGaussianParams{"CL2Sigma", kCL2Sigma}, + PlcGaussianParams{"CL1Sigma", kCL1Sigma}, + PlcGaussianParams{"CLNearZero", 0.00001})), // boundary case CL -> 0 + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// PROFILE LIKELIHOOD CALCULATOR - LIKELIHOOD INTERVAL - POISSON DISTRIBUTION +// +// Test the 68% likelihood interval computed by the profile likelihood +// calculator on a Poisson distribution, from only one observed value. +// Reference values are computed via analytic methods: solve equation +// 2*[ln(LL(xMax)) - ln(LL(x))] = 1. +// +// ModelConfig (implicit) : +// Observable -> x +// Parameter of Interest -> mean +// +// Original tests 6-10 (TestProfileLikelihoodCalculator2), with the observed +// value probing the boundaries of the [0,1000] range. +// +/////////////////////////////////////////////////////////////////////////////// + +struct PlcPoissonParams { + std::string name; + int obsValue; +}; + +class PlcPoissonInterval : public StressRooStatsParamTest {}; + +TEST_P(PlcPoissonInterval, CompareWithAnalyticInterval) +{ + const int obsValue = param().obsValue; + + // Set a 68% confidence level for the interval + const double confidenceLevel = kCL1Sigma; + + // Create Poisson model and dataset + RooWorkspace ws{"w"}; + ws.factory(TString::Format("Poisson::poiss(x[%d,0,1000], mean[0,1000])", obsValue).Data()); + RooDataSet data{"data", "data", *ws.var("x")}; + data.add(*ws.var("x")); + + // Calculate likelihood interval using the ProfileLikelihoodCalculator + ProfileLikelihoodCalculator plc{data, *ws.pdf("poiss"), *ws.var("mean")}; + plc.SetConfidenceLevel(confidenceLevel); + std::unique_ptr interval{plc.GetInterval()}; + + // Reference limits are the solutions of 2*[ln(LL(xMax)) - ln(LL(x))] = 1, + // where xMax is the point of maximum likelihood. For the special case of + // the Poisson distribution with N = 1, xMax = obsValue. + TString llRatioExpression = + TString::Format("2*(x-%d*log(x)-%d+%d*log(%d))", obsValue, obsValue, obsValue, obsValue); + // Special case obsValue = 0 because log(0) is not computable, the limit of + // n * log(n), n->0 must be taken + if (obsValue == 0) + llRatioExpression = "2*x"; + TF1 llRatio{"llRatio", llRatioExpression, 1e-100, double(obsValue)}; // lowerLimit < obsValue + + // For obsValue = 0 there is no analytic lower limit (the reference value in + // the original suite was NaN, which made its comparison pass trivially) + if (obsValue != 0) { + EXPECT_NEAR(interval->LowerLimit(*ws.var("mean")), llRatio.GetX(1), kVTol); } + llRatio.SetRange(obsValue, 1000); // upperLimit > obsValue + EXPECT_NEAR(interval->UpperLimit(*ws.var("mean")), llRatio.GetX(1), kVTol); +} - delete gBenchmark; - gBenchmark = nullptr; +INSTANTIATE_TEST_SUITE_P(RooStats, PlcPoissonInterval, + ::testing::Combine(::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values( + PlcPoissonParams{"Obs0", 0}, // boundary Poisson value (0) + PlcPoissonParams{"Obs1", 1}, PlcPoissonParams{"Obs5", 5}, + PlcPoissonParams{"Obs100", 100}, + PlcPoissonParams{"Obs800", 800})), // boundary Poisson value + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// PROFILE LIKELIHOOD CALCULATOR - LIKELIHOOD INTERVAL - POISSON PRODUCT MODEL +// +// Test the 68% likelihood interval computed by the ProfileLikelihoodCalculator +// on a complex model. Reference values and test values are both computed with +// the ProfileLikelihoodCalculator. As such, this test can only confirm if the +// ProfileLikelihoodCalculator has the same behaviour across different computer +// platforms or RooStats revisions. +// +// ModelConfig (explicit) : Poisson Product Model +// built in stressRooStats_models.h +// +// Original tests 11-13 (TestProfileLikelihoodCalculator3). +// +/////////////////////////////////////////////////////////////////////////////// + +struct PoissonProductParams { + std::string name; + int obsValueX; // observed value "x" when measuring sig + bkg1 + int obsValueY; // observed value "y" when measuring 2*sig*1.2^beta + bkg2 + double confidenceLevel; + double refLowerLimit; + double refUpperLimit; +}; + +class PlcPoissonProductInterval : public StressRooStatsParamTest {}; + +TEST_P(PlcPoissonProductInterval, RegressionInterval) +{ + // Create workspace and model with the observed values in the data set + RooWorkspace ws{"w"}; + ModelConfig *model = setupPoissonProductModel(ws, param().obsValueX, param().obsValueY); + ASSERT_NE(model, nullptr); + + // build likelihood interval with ProfileLikelihoodCalculator + ProfileLikelihoodCalculator plc{*ws.data("data"), *model}; + plc.SetConfidenceLevel(param().confidenceLevel); + std::unique_ptr interval{plc.GetInterval()}; + + EXPECT_NEAR(interval->LowerLimit(*ws.var("sig")), param().refLowerLimit, kVTol); + EXPECT_NEAR(interval->UpperLimit(*ws.var("sig")), param().refUpperLimit, kVTol); +} - // Some of the object are multiple times in the list, let's make sure they - // are not deleted twice. - // The addition of memDir to the list of Cleanups is not needed if it already - // there, for example if memDir is gROOT. - bool needCleanupAdd = nullptr == gROOT->GetListOfCleanups()->FindObject(memDir->GetList()); - if (needCleanupAdd) - gROOT->GetListOfCleanups()->Add(memDir->GetList()); +INSTANTIATE_TEST_SUITE_P( + RooStats, PlcPoissonProductInterval, + ::testing::Combine( + ::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values( + PoissonProductParams{"Obs10_30_CL1Sigma", 10, 30, kCL1Sigma, 7.1428311212946145, 12.216799640536237}, + PoissonProductParams{"Obs20_25_CL1Sigma", 20, 25, kCL1Sigma, 8.7079993530346087, 14.911357163771072}, + PoissonProductParams{"Obs15_20_CL2Sigma", 15, 20, kCL2Sigma, 3.4494551349741052, 14.067570959270299})), + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// PROFILE LIKELIHOOD CALCULATOR HYPOTHESIS TEST - ON / OFF MODEL +// +// Perform a hypothesis test using the ProfileLikelihoodCalculator on the +// on/off model. The reference values are taken from the paper: "Evaluation +// of three methods for calculating statistical significance when incorporating +// a systematic uncertainty into a test of the background-only hypothesis for +// a Poisson process" by Robert D. Cousins, James T. Linnemann, Jordan Tucker. +// +// ModelConfig (explicit) : Poisson On / Off Model +// built in stressRooStats_models.h +// +// Original test 14 (TestProfileLikelihoodCalculator4). +// +/////////////////////////////////////////////////////////////////////////////// + +class PlcOnOffHypoTest : public StressRooStatsBackendTest {}; + +TEST_P(PlcOnOffHypoTest, CompareSignificanceWithPaperValues) +{ + // A larger tolerance is needed since the values in the Cousins paper are + // given with 1e-2 precision + const double tol = 1e-2; + + // For testing purposes, we consider three special cases for which the + // values are known from the Cousins et al. paper mentioned above. The + // inputs for each of these cases are (using the notations from the paper): + // n_on, n_off, tau and Z_PL. + const int numberTestSets = 3; + const int numberOnEvents[numberTestSets] = {4, 50, 67}; + const int numberOffEvents[numberTestSets] = {5, 55, 15}; + const double tau[numberTestSets] = {5.0, 2.0, 0.5}; + const double significance[numberTestSets] = {1.95, 3.02, 3.04}; + + for (int i = 0; i < numberTestSets; ++i) { + + // build workspace and model + RooWorkspace ws{"w"}; + buildOnOffModel(ws); + auto sbModel = dynamic_cast(ws.obj("S+B")); + auto bModel = dynamic_cast(ws.obj("B")); + ASSERT_NE(sbModel, nullptr); + ASSERT_NE(bModel, nullptr); + + // add observable values to data set + ws.var("n_on")->setVal(numberOnEvents[i]); + ws.var("n_off")->setVal(numberOffEvents[i]); + ws.var("tau")->setVal(tau[i]); + ws.var("tau")->setConstant(); + ws.data("data")->add(*sbModel->GetObservables()); + + // set snapshots + ws.var("sig")->setVal(numberOnEvents[i] - numberOffEvents[i] / tau[i]); + sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); + ws.var("sig")->setVal(0); + bModel->SetSnapshot(*bModel->GetParametersOfInterest()); + + // has as initial value a non-zero value for sig (i.e start with the S+B value) + sbModel->LoadSnapshot(); + + // get significance using the ProfileLikelihoodCalculator + ProfileLikelihoodCalculator plc{*ws.data("data"), *sbModel}; + plc.SetNullParameters(*bModel->GetSnapshot()); + + std::unique_ptr result{plc.GetHypoTest()}; + EXPECT_NEAR(result->Significance(), significance[i], tol) + << "for n_on = " << numberOnEvents[i] << ", n_off = " << numberOffEvents[i] << ", tau = " << tau[i]; + } +} - memDir->GetList()->Delete("slow"); +INSTANTIATE_TEST_SUITE_P(RooStats, PlcOnOffHypoTest, ::testing::Values(ROOFIT_EVAL_BACKENDS), backendName); + +/////////////////////////////////////////////////////////////////////////////// +// +// PART TWO: BAYESIAN CALCULATOR UNIT TESTS +// +/////////////////////////////////////////////////////////////////////////////// + +/////////////////////////////////////////////////////////////////////////////// +// +// BAYESIAN CENTRAL INTERVAL - SIMPLE POISSON MODEL +// +// Test the Bayesian central interval computed by the BayesianCalculator on a +// Poisson distribution, using different priors. The parameter of interest is +// the mean of the Poisson distribution, and there are no nuisance parameters. +// The priors used are: +// 1. constant / uniform +// 2. inverse of the mean +// 3. square root of the inverse of the mean +// 4. gamma distribution +// The posterior distribution is easily obtained analytically for these cases. +// Therefore, the reference interval limits are computed analytically. +// +// ModelConfig (implicit) : +// Observable -> x +// Parameter of Interest -> mean +// +// Original tests 15-18 (TestBayesianCalculator1). +// +/////////////////////////////////////////////////////////////////////////////// + +namespace { +double priorInvFunc(double mean) +{ + return 1.0 / mean; +} +double priorInvSqrtFunc(double mean) +{ + return 1.0 / std::sqrt(mean); +} +} // namespace + +struct BcPoissonParams { + std::string name; + int obsValue; +}; - if (needCleanupAdd) - gROOT->GetListOfCleanups()->Remove(memDir->GetList()); +class BcCentralInterval : public StressRooStatsParamTest {}; - return nFailed; +TEST_P(BcCentralInterval, CompareWithAnalyticInterval) +{ + const int obsValue = param().obsValue; + + // Set the confidence level for a 68.3% CL central interval + const double confidenceLevel = kCL1Sigma; + const double gammaShape = 2; // shape of the gamma distribution prior (gamma = alpha) + const double gammaRate = 1; // rate = 1/scale of the gamma distribution prior (beta = 1/theta) + const int numberScans = 10000; // tested to be sufficient for the scan of the Bayesian posterior + + // Create Poisson model + RooWorkspace ws{"w"}; + ws.factory("Poisson::poiss(x[0,100], mean[1e-6,100])"); + + // create prior pdfs + ws.factory("Uniform::prior(mean)"); + ws.import(RooCFunction1PdfBinding("priorInv", "priorInv", &priorInvFunc, *ws.var("mean"))); + ws.import( + RooCFunction1PdfBinding("priorInvSqrt", "priorInvSqrt", &priorInvSqrtFunc, *ws.var("mean"))); + ws.factory(TString::Format("Gamma::priorGamma(mean, %lf, %lf, 0)", gammaShape, gammaRate).Data()); + + // build argument sets and data set + ws.defineSet("obs", "x"); + ws.defineSet("poi", "mean"); + ws.var("x")->setVal(obsValue); + RooDataSet data{"data", "data", *ws.set("obs")}; + data.add(*ws.set("obs")); + + // Compute the interval with the BayesianCalculator for the given prior and + // compare with the analytically computed reference limits: the posterior + // for a Poisson model with the priors used here is a gamma distribution. + auto testPrior = [&](const char *priorName, double refLowerLimit, double refUpperLimit) { + BayesianCalculator bc{data, *ws.pdf("poiss"), *ws.set("poi"), *ws.pdf(priorName), nullptr}; + bc.SetConfidenceLevel(confidenceLevel); + bc.SetScanOfPosterior(numberScans); + std::unique_ptr interval{bc.GetInterval()}; + EXPECT_NEAR(interval->LowerLimit(), refLowerLimit, kVTol) << "lower limit for prior " << priorName; + EXPECT_NEAR(interval->UpperLimit(), refUpperLimit, kVTol) << "upper limit for prior " << priorName; + }; + + const double testSize = (1.0 - confidenceLevel) / 2; + + // Uniform prior on mean + testPrior("prior", gamma_quantile(testSize, obsValue + 1, 1), // integrate to 16% + gamma_quantile_c(testSize, obsValue + 1, 1)); // integrate to 84% + // Inverse of mean prior + testPrior("priorInv", gamma_quantile(testSize, obsValue, 1), gamma_quantile_c(testSize, obsValue, 1)); + // Square root of inverse of mean prior + testPrior("priorInvSqrt", gamma_quantile(testSize, obsValue + 0.5, 1), + gamma_quantile_c(testSize, obsValue + 0.5, 1)); + // Gamma distribution prior + testPrior("priorGamma", gamma_quantile(testSize, obsValue + gammaShape, 1.0 / (1 + gammaRate)), + gamma_quantile_c(testSize, obsValue + gammaShape, 1.0 / (1 + gammaRate))); } -//_____________________________batch only_____________________ -#ifndef __CLING__ +INSTANTIATE_TEST_SUITE_P(RooStats, BcCentralInterval, + ::testing::Combine(::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values(BcPoissonParams{"Obs1", 1}, BcPoissonParams{"Obs3", 3}, + BcPoissonParams{"Obs10", 10}, + BcPoissonParams{"Obs50", 50})), + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// BAYESIAN SHORTEST INTERVAL - SIMPLE POISSON MODEL +// +// Test the Bayesian shortest interval computed by the BayesianCalculator on a +// Poisson distribution, using different priors. The parameter of interest is +// the mean of the Poisson distribution, and there are no nuisance parameters. +// The priors used are: +// 1. constant / uniform +// 2. inverse of the mean +// The reference interval limits are taken from the paper: "Why isn't every +// physicist a Bayesian?" by Robert D. Cousins. +// +// Original test 19 (TestBayesianCalculator2). +// +/////////////////////////////////////////////////////////////////////////////// + +class BcShortestInterval : public StressRooStatsBackendTest {}; + +TEST_P(BcShortestInterval, CompareWithPaperValues) +{ + // the reference values in the paper have a precision of only two decimal + // points, so we increase the value tolerance accordingly + const double tol = 1e-2; + + // Put the confidence level so that we obtain a 68% confidence interval + const double confidenceLevel = kCL1Sigma; + const int obsValue = 3; // observed experiment value + const int numberScans = 100000; // sufficient number of scans + + // Create Poisson model + RooWorkspace ws{"w"}; + ws.factory("Poisson::poiss(x[0,100], mean[1e-6,100])"); + ws.factory("Uniform::prior(mean)"); + ws.factory("EXPR::priorInv('1/mean', mean)"); + + // build argument sets and data set + ws.defineSet("poi", "mean"); + ws.defineSet("obs", "x"); + ws.var("x")->setVal(obsValue); + RooDataSet data{"data", "data", *ws.set("obs")}; + data.add(*ws.set("obs")); + + auto testPrior = [&](const char *priorName, double refLowerLimit, double refUpperLimit) { + BayesianCalculator bc{data, *ws.pdf("poiss"), *ws.set("poi"), *ws.pdf(priorName), nullptr}; + bc.SetConfidenceLevel(confidenceLevel); + bc.SetShortestInterval(); + bc.SetScanOfPosterior(numberScans); + std::unique_ptr interval{bc.GetInterval()}; + EXPECT_NEAR(interval->LowerLimit(), refLowerLimit, tol) << "lower limit for prior " << priorName; + EXPECT_NEAR(interval->UpperLimit(), refUpperLimit, tol) << "upper limit for prior " << priorName; + }; + + // Uniform prior on mean + testPrior("prior", 1.55, 5.15); + // Inverse of mean prior + testPrior("priorInv", 0.86, 3.85); +} -int main(int argc, const char *argv[]) +INSTANTIATE_TEST_SUITE_P(RooStats, BcShortestInterval, ::testing::Values(ROOFIT_EVAL_BACKENDS), backendName); + +/////////////////////////////////////////////////////////////////////////////// +// +// BAYESIAN CENTRAL INTERVAL - POISSON PRODUCT MODEL +// +// Test the validity of the central interval computed by the BayesianCalculator +// on a complex Poisson model distribution. Reference values and test values +// are both computed with the BayesianCalculator. As such, this test can only +// confirm if the BayesianCalculator has the same behaviour across different +// computing platforms or RooStats revisions. A uniform prior PDF is used for +// the parameter of interest ("sig"). +// +// ModelConfig (explicit) : Poisson Product Model +// built in stressRooStats_models.h +// +// Original tests 20-22 (TestBayesianCalculator3). +// +/////////////////////////////////////////////////////////////////////////////// + +class BcPoissonProductInterval : public StressRooStatsParamTest {}; + +TEST_P(BcPoissonProductInterval, RegressionInterval) { - bool doWrite = false; - int verbose = 0; - bool allTests = false; - bool oneTest = false; - int testNumber = 0; - bool dryRun = false; - bool doDump = false; - bool doTreeStore = false; - auto backend = RooFit::EvalBackend(RooFit::EvalBackend::Value::Legacy); - - // string refFileName = "http://root.cern/files/stressRooStats_v534_ref.root" ; - string refFileName = "stressRooStats_ref.root"; - string minimizerName = "Minuit"; - - // Parse command line arguments - for (int i = 1; i < argc; i++) { - string arg = argv[i]; - - if (arg == "-b") { - std::string mode = argv[++i]; - backend = RooFit::EvalBackend(mode); - std::cout << "stressRooStats: NLL evaluation backend set to " << mode << std::endl; - } else if (arg == "-f") { - std::cout << "stressRooStats: using reference file " << argv[i + 1] << std::endl; - refFileName = argv[++i]; - } else if (arg == "-w") { - std::cout << "stressRooStats: running in writing mode to update reference file" << std::endl; - doWrite = true; - } else if (arg == "-mc") { - std::cout << "stressRooStats: running in memcheck mode, no regression tests are performed" << std::endl; - dryRun = true; - } else if (arg == "-min" || arg == "-minim") { - std::cout << "stressRooStats: running using minimizer " << argv[i + 1] << std::endl; - minimizerName = argv[++i]; - } else if (arg == "-ts") { - std::cout << "stressRooStats: setting tree-based storage for datasets" << std::endl; - doTreeStore = true; - } else if (arg == "-v") { - std::cout << "stressRooStats: running in verbose mode" << std::endl; - verbose = 1; - } else if (arg == "-vv") { - std::cout << "stressRooStats: running in very verbose mode" << std::endl; - verbose = 2; - } else if (arg == "-vvv") { - std::cout << "stressRooStats: running in very very verbose mode" << std::endl; - verbose = 3; - } else if (arg == "-a") { - std::cout << "stressRooStats: deploying full suite of tests" << std::endl; - allTests = true; - } else if (arg == "-n") { - std::cout << "stressRooStats: running single test" << std::endl; - oneTest = true; - testNumber = atoi(argv[++i]); - } else if (arg == "-d") { - std::cout << "stressRooStats: setting gDebug to " << argv[i + 1] << std::endl; - gDebug = atoi(argv[++i]); - } else if (arg == "-c") { - std::cout << "stressRooStats: dumping comparison file for failed tests " << std::endl; - doDump = true; - } else if (arg == "-h" || arg == "--help") { - std::cout << R"(usage: stressRooStats [ options ] - - -b : Perform every fit in the tests with the EvalBackend() command argument, where is a string - -f : use given reference file instead of default ("stressRooStats_ref.root") - -w : write reference file, instead of reading file and running comparison tests - -n N : only run test with sequential number N - -a : run full suite of tests (default is basic suite); this overrides the -n single test option - -c : dump file stressRooStats_DEBUG.root to which results of both current result and reference for each failed test are written - -mc : memory check mode, no regression test are performed. Set this flag when running with valgrind - -min : minimizer name (default is Minuit, not Minuit2) - -vs : use vector-based storage for all datasets (default is tree-based storage) - -v/-vv : set verbose mode (show result of each regression test) or very verbose mode (show all roofit output as well) - -d N : set ROOT gDebug flag to N -)"; - return 0; - } + const int numberScans = 10; // sufficient number of scans + + // Create workspace and model with the observed values in the data set + RooWorkspace ws{"w"}; + ModelConfig *model = setupPoissonProductModel(ws, param().obsValueX, param().obsValueY); + ASSERT_NE(model, nullptr); + + // Create BayesianCalculator + BayesianCalculator bc{*ws.data("data"), *model}; + bc.SetConfidenceLevel(param().confidenceLevel); + bc.SetScanOfPosterior(numberScans); + + // Obtain confidence interval by scanning the posterior function in the + // given number of points + std::unique_ptr interval{bc.GetInterval()}; + EXPECT_NEAR(interval->LowerLimit(), param().refLowerLimit, kVTol); + EXPECT_NEAR(interval->UpperLimit(), param().refUpperLimit, kVTol); +} + +INSTANTIATE_TEST_SUITE_P( + RooStats, BcPoissonProductInterval, + ::testing::Combine( + ::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values( + PoissonProductParams{"Obs10_30_CL1Sigma", 10, 30, kCL1Sigma, 7.1665080051981258, 12.312785237149114}, + PoissonProductParams{"Obs20_25_CL1Sigma", 20, 25, kCL1Sigma, 8.7831170668504654, 14.874961130060584}, + PoissonProductParams{"Obs15_20_CL2Sigma", 15, 20, kCL2Sigma, 3.4603352565856857, 14.186182799724543})), + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// PART THREE: MARKOV CHAIN MONTE CARLO CALCULATOR UNIT TESTS +// +/////////////////////////////////////////////////////////////////////////////// + +/////////////////////////////////////////////////////////////////////////////// +// +// MCMC INTERVAL CALCULATOR - POISSON PRODUCT MODEL +// +// Test the validity of the confidence interval computed by the MCMCCalculator +// on a complex Poisson model distribution. Reference values and test values +// are both computed with the MCMCCalculator. As such, this test can only +// confirm if the MCMCCalculator has the same behaviour across different +// computing platforms or RooStats revisions. +// +// ModelConfig (explicit) : Poisson Product Model +// built in stressRooStats_models.h +// +// Original tests 23-25 (TestMCMCCalculator). +// +/////////////////////////////////////////////////////////////////////////////// + +class McmcPoissonProductInterval : public StressRooStatsParamTest {}; + +TEST_P(McmcPoissonProductInterval, RegressionInterval) +{ + // Create workspace and model with the observed values in the data set + RooWorkspace ws{"w"}; + ModelConfig *model = setupPoissonProductModel(ws, param().obsValueX, param().obsValueY); + ASSERT_NE(model, nullptr); + + // create and configure MCMC calculator + SequentialProposal sp{0.1}; + MCMCCalculator mcmcc{*ws.data("data"), *model}; + mcmcc.SetProposalFunction(sp); + mcmcc.SetNumIters(100000); // Metropolis-Hastings algorithm iterations + mcmcc.SetNumBurnInSteps(50); // first 50 steps to be ignored as burn-in + mcmcc.SetConfidenceLevel(param().confidenceLevel); + + // calculate the confidence interval + std::unique_ptr interval{mcmcc.GetInterval()}; + EXPECT_NEAR(interval->LowerLimit(*ws.var("sig")), param().refLowerLimit, kVTol); + EXPECT_NEAR(interval->UpperLimit(*ws.var("sig")), param().refUpperLimit, kVTol); +} + +INSTANTIATE_TEST_SUITE_P( + RooStats, McmcPoissonProductInterval, + ::testing::Combine(::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values(PoissonProductParams{"Obs10_30_CL1Sigma", 10, 30, kCL1Sigma, 7.0, 11.4}, + PoissonProductParams{"Obs20_25_CL1Sigma", 20, 25, kCL1Sigma, 9.0, + 14.600000000000001}, + PoissonProductParams{"Obs15_20_CL2Sigma", 15, 20, kCL2Sigma, + 3.4000000000000004, 13.800000000000001})), + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// PART FOUR: HYPOTHESIS TEST CALCULATOR UNIT TESTS +// +/////////////////////////////////////////////////////////////////////////////// + +/////////////////////////////////////////////////////////////////////////////// +// +// ZBI - ON / OFF MODEL +// +// Evaluate the functionality of the top level function +// NumberCountingUtils::BinomialWithTauObsZ, which computes the significance +// of a hypothesis test via a frequentist solution. This significance, called +// ZBi, is detailed in the article "Evaluation of three methods for calculating +// statistical significance when incorporating a systematic uncertainty into a +// test of the background-only hypothesis for a Poisson process" by Robert D. +// Cousins, James T. Linnemann, Jordan Tucker. The reference values are taken +// from the paper. +// +// This computation involves no likelihood fits, so it is not parameterized +// over the evaluation backends. +// +// Original test 26 (TestZBi). +// +/////////////////////////////////////////////////////////////////////////////// + +TEST(ZBiSignificance, CompareWithPaperValues) +{ + setUpStressRooStatsTest(); + + // A larger tolerance is needed since the values in the Cousins paper are + // given with 1e-2 precision + const double tol = 1e-2; + + // For testing purposes, we consider four special cases for which the values + // are known from the Cousins et al. paper mentioned above. The inputs for + // each of these cases are (using the notations from the paper): n_on, n_off + // and tau. + const int numberTestSets = 4; + const int numberOnEvents[numberTestSets] = {4, 50, 67, 200}; + const int numberOffEvents[numberTestSets] = {5, 55, 15, 10}; + const double tau[numberTestSets] = {5.0, 2.0, 0.5, 0.1}; + const double significance[numberTestSets] = {1.66, 2.93, 2.89, 2.2}; + + for (int i = 0; i < numberTestSets; ++i) { + EXPECT_NEAR(NumberCountingUtils::BinomialWithTauObsZ(numberOnEvents[i], numberOffEvents[i], tau[i]), + significance[i], tol) + << "for n_on = " << numberOnEvents[i] << ", n_off = " << numberOffEvents[i] << ", tau = " << tau[i]; } - // if (doWrite && refFileName.find("http:") == 0) { + tearDownStressRooStatsTest(); +} - // // Locate file name part in URL and update refFileName accordingly - // char* buf = new char[refFileName.size() + 1]; - // strcpy(buf, refFileName.c_str()); - // char *ptr = strrchr(buf, '/'); - // if (!ptr) ptr = strrchr(buf, ':'); - // refFileName = ptr + 1; - // delete[] buf; +/////////////////////////////////////////////////////////////////////////////// +// +// ASYMPTOTIC CALCULATOR VS PROFILE LIKELIHOOD CALCULATOR HYPOTHESIS TEST +// +// This test evaluates the functionality of the AsymptoticCalculator by +// comparing the significance given from a hypothesis test on the on/off model +// with the significance given by the ProfileLikelihoodCalculator. If working +// properly, the two methods should yield identical results. On top of the +// direct comparison of the two methods, both significances are also compared +// with the frozen reference values of the original suite (which were produced +// with the ProfileLikelihoodCalculator). +// +// ModelConfig (explicit) : Poisson On / Off Model +// built in stressRooStats_models.h +// +// Original tests 27-31 (TestHypoTestCalculator1). +// +/////////////////////////////////////////////////////////////////////////////// + +struct OnOffSignificanceParams { + std::string name; + int obsValueOn; // observed value "n_on" of sig + bkg + int obsValueOff; // observed value "n_off" of tau * bkg + double tau; // parameter of the model (constant with regard to integration) + double refSignificance; +}; + +class AsymptoticVsPlcSignificance : public StressRooStatsParamTest {}; + +TEST_P(AsymptoticVsPlcSignificance, CompareSignificances) +{ + // build workspace and model, add observable values to the data set and fix + // other parameters, then make the S+B and B snapshots + auto setupWorkspace = [&](RooWorkspace &ws) -> std::pair { + buildOnOffModel(ws); + auto sbModel = dynamic_cast(ws.obj("S+B")); + auto bModel = dynamic_cast(ws.obj("B")); + if (!sbModel || !bModel) + return {nullptr, nullptr}; + + ws.var("n_on")->setVal(param().obsValueOn); + ws.var("n_off")->setVal(param().obsValueOff); + ws.var("tau")->setVal(param().tau); + ws.var("tau")->setConstant(); + ws.data("data")->add(*sbModel->GetObservables()); + ws.var("bkg")->setVal(param().obsValueOff / param().tau); + + ws.var("sig")->setVal(param().obsValueOn - param().obsValueOff / param().tau); + sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); + ws.var("sig")->setVal(0.0); + bModel->SetSnapshot(*bModel->GetParametersOfInterest()); + + return std::make_pair(sbModel, bModel); + }; + + // Hypothesis test with the ProfileLikelihoodCalculator + double significancePlc = 0.0; + { + RooWorkspace ws{"w"}; + auto [sbModel, bModel] = setupWorkspace(ws); + ASSERT_NE(sbModel, nullptr); + ASSERT_NE(bModel, nullptr); + + ProfileLikelihoodCalculator plc{*ws.data("data"), *sbModel}; + plc.SetNullParameters(*bModel->GetSnapshot()); + plc.SetAlternateParameters(*sbModel->GetSnapshot()); + std::unique_ptr result{plc.GetHypoTest()}; + significancePlc = result->Significance(); + } - // std::cout << "stressRooStats: WARNING running in write mode, but reference file is web file, writing local file - // instead: " - // << refFileName << std::endl; - // } + // Hypothesis test with the AsymptoticCalculator + double significanceAc = 0.0; + { + RooWorkspace ws{"w"}; + auto [sbModel, bModel] = setupWorkspace(ws); + ASSERT_NE(sbModel, nullptr); + ASSERT_NE(bModel, nullptr); + + AsymptoticCalculator atc{*ws.data("data"), *sbModel, *bModel}; + atc.SetOneSidedDiscovery(true); + std::unique_ptr result{atc.GetHypoTest()}; + significanceAc = result->Significance(); + } - // set minimizer - ROOT::Math::MinimizerOptions::SetDefaultMinimizer(minimizerName.c_str()); + EXPECT_NEAR(significancePlc, param().refSignificance, kVTol); + EXPECT_NEAR(significanceAc, param().refSignificance, kVTol); +} - // set default NLL backend - RooFit::EvalBackend::defaultValue() = backend.value(); +INSTANTIATE_TEST_SUITE_P( + RooStats, AsymptoticVsPlcSignificance, + ::testing::Combine( + ::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values(OnOffSignificanceParams{"Obs150_100_Tau1", 150, 100, 1.0, 3.1729725711197787}, + OnOffSignificanceParams{"Obs200_100_Tau1", 200, 100, 1.0, 5.8292198861163191}, + OnOffSignificanceParams{"Obs105_100_Tau1", 105, 100, 1.0, 0.34923144876992013}, + OnOffSignificanceParams{"Obs150_10_Tau0p1", 150, 10, 0.1, 1.3181915918803733}, + OnOffSignificanceParams{"Obs150_400_Tau4", 150, 400, 4.0, 4.0985771507906898})), + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// HYPOTHESIS TEST CALCULATOR TEST - SIMULTANEOUS PDF MODEL +// +// This test evaluates the functionality of the HypoTestCalculator by +// calculating the significance of the signal on a simple Simultaneous Pdf +// model with two channels. Reference values and test values are both computed +// with the HypoTestCalculator. As such, this test can only confirm if the +// HypoTestCalculator has the same behaviour across different computing +// platforms or RooStats revisions. +// +// ModelConfig (explicit) : Simultaneous Model +// built in stressRooStats_models.h +// +// Original tests 32-36 (TestHypoTestCalculator2). +// +/////////////////////////////////////////////////////////////////////////////// + +struct HypoTestCalculatorParams { + std::string name; + ECalculatorType calculatorType; + ETestStatType testStatType; + double refSignificance; +}; + +class HypoTestCalculatorSignificance : public StressRooStatsParamTest {}; + +TEST_P(HypoTestCalculatorSignificance, RegressionSignificance) +{ + // Build workspace and models + RooWorkspace ws{"w"}; + buildSimultaneousModel(&ws); + auto sbModel = dynamic_cast(ws.obj("S+B")); + auto bModel = dynamic_cast(ws.obj("B")); + ASSERT_NE(sbModel, nullptr); + ASSERT_NE(bModel, nullptr); + + // set snapshots + sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); // value set in model + ws.var("sig")->setVal(0); + bModel->SetSnapshot(*bModel->GetParametersOfInterest()); + + // the test statistic is declared before the calculator because the + // calculator's sampler will hold a raw pointer to it + std::unique_ptr testStat{buildTestStatistic(param().testStatType, *bModel, *sbModel)}; + std::unique_ptr calc{ + buildHypoTestCalculator(param().calculatorType, *ws.data("data"), *bModel, *sbModel, 500, 50)}; + if (param().calculatorType == kAsymptotic) { + static_cast(*calc).SetOneSidedDiscovery(true); + } + + // ToyMCSampler configuration + auto tmcs = static_cast(calc->GetTestStatSampler()); + tmcs->SetTestStatistic(testStat.get()); + tmcs->SetUseMultiGen(true); // speedup - gBenchmark = new TBenchmark(); - return stressRooStats(refFileName.c_str(), doWrite, verbose, allTests, oneTest, testNumber, dryRun, doDump, - doTreeStore); + std::unique_ptr result{calc->GetHypoTest()}; + EXPECT_NEAR(result->Significance(), param().refSignificance, kVTol); } -//////////////////////////////////////////////////////////////////////////////// +INSTANTIATE_TEST_SUITE_P( + RooStats, HypoTestCalculatorSignificance, + ::testing::Combine( + ::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values(HypoTestCalculatorParams{"Asymptotic_ProfileLROneSidedDiscovery", kAsymptotic, + kProfileLROneSidedDiscovery, 2.2499193885237001}, + HypoTestCalculatorParams{"Frequentist_SimpleLR", kFrequentist, kSimpleLR, 2.3263478740408408}, + HypoTestCalculatorParams{"Frequentist_RatioLR", kFrequentist, kRatioLR, 2.4089155458154612}, + HypoTestCalculatorParams{"Frequentist_ProfileLROneSidedDiscovery", kFrequentist, + kProfileLROneSidedDiscovery, 2.5121443279304616}, + HypoTestCalculatorParams{"Hybrid_ProfileLROneSidedDiscovery", kHybrid, + kProfileLROneSidedDiscovery, 2.2571292444862254})), + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// PART FIVE: HYPOTHESIS TEST INVERTER UNIT TESTS +// +/////////////////////////////////////////////////////////////////////////////// + +/////////////////////////////////////////////////////////////////////////////// +// +// HYPOTESTINVERTER INTERVAL - POISSON PRODUCT MODEL +// +// Test the validity of the confidence interval computed by the +// HypoTestInverter on a complex Poisson model distribution. Reference values +// and test values are both computed with the HypoTestInverter. As such, this +// test can only confirm if the HypoTestInverter has the same behaviour across +// different computing platforms or RooStats revisions. +// +// ModelConfig (explicit) : Poisson Product Model +// built in stressRooStats_models.h +// +// Original tests 37-43 (TestHypoTestInverter1). +// +/////////////////////////////////////////////////////////////////////////////// + +struct HypoTestInverterParams { + std::string name; + ECalculatorType calculatorType; + ETestStatType testStatType; + int obsValueX; // observed value "x" when measuring sig + bkg1 + int obsValueY; // observed value "y" when measuring 2*sig*1.2^beta + bkg2 + double refLowerLimit; + double refUpperLimit; +}; + +class HypoTestInverterInterval : public StressRooStatsParamTest {}; + +TEST_P(HypoTestInverterInterval, RegressionInterval) +{ + const double confidenceLevel = kCL1Sigma; + + // larger value test tolerance especially when using toys (difference of + // <~ 0.1 observed between using Minuit or Minuit2) + const double tol = (param().calculatorType == kAsymptotic) ? 0.01 : 0.1; + + // Create workspace and model with the observed values in the data set + RooWorkspace ws{"w"}; + ModelConfig *sbModel = setupPoissonProductModel(ws, param().obsValueX, param().obsValueY); + auto bModel = dynamic_cast(ws.obj("B")); + ASSERT_NE(sbModel, nullptr); + ASSERT_NE(bModel, nullptr); + + // set snapshots + ws.var("sig")->setVal(param().obsValueX - ws.var("bkg1")->getValV()); + sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); + ws.var("sig")->setVal(0); + bModel->SetSnapshot(*bModel->GetParametersOfInterest()); + + // build and configure HypoTestInverter (the test statistic is declared + // first because the calculator will hold a raw pointer to it) + std::unique_ptr testStat{buildTestStatistic(param().testStatType, *sbModel, *bModel)}; + std::unique_ptr calc{ + buildHypoTestCalculator(param().calculatorType, *ws.data("data"), *sbModel, *bModel, 100, 1)}; + HypoTestInverter hti{*calc, nullptr, 1.0 - confidenceLevel}; + hti.SetTestStatistic(*testStat); + + int nscanPoints = 10; + if (param().calculatorType == kAsymptotic) { + static_cast(*calc).SetTwoSided(); + nscanPoints = 40; + } + + hti.SetFixedScan(nscanPoints, ws.var("sig")->getMin(), ws.var("sig")->getMax()); // significant speedup + + // ToyMCSampler configuration + auto tmcs = static_cast(hti.GetHypoTestCalculator()->GetTestStatSampler()); + tmcs->SetNEventsPerToy(1); // needed because we don't have an extended pdf + tmcs->SetUseMultiGen(true); // speedup -int stressRooStats() + std::unique_ptr interval{hti.GetInterval()}; + EXPECT_NEAR(interval->LowerLimit(), param().refLowerLimit, tol); + EXPECT_NEAR(interval->UpperLimit(), param().refUpperLimit, tol); +} + +INSTANTIATE_TEST_SUITE_P( + RooStats, HypoTestInverterInterval, + ::testing::Combine(::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values(HypoTestInverterParams{"Asymptotic_ProfileLR_Obs10_30", kAsymptotic, + kProfileLR, 10, 30, 7.1386012828228216, + 12.223793911203613}, + HypoTestInverterParams{"Asymptotic_ProfileLR_Obs20_25", kAsymptotic, + kProfileLR, 20, 25, 8.7069484829792732, + 14.91407242494005}, + HypoTestInverterParams{"Asymptotic_ProfileLR_Obs15_20", kAsymptotic, + kProfileLR, 15, 20, 5.750533033549976, + 10.948013504988015}, + HypoTestInverterParams{"Frequentist_ProfileLR_Obs10_30", kFrequentist, + kProfileLR, 10, 30, 7.2115486004250906, + 12.341346144450787}, + HypoTestInverterParams{"Frequentist_ProfileLR_Obs20_25", kFrequentist, + kProfileLR, 20, 25, 8.6294497270104724, + 15.078331964885663}, + HypoTestInverterParams{"Frequentist_ProfileLR_Obs15_20", kFrequentist, + kProfileLR, 15, 20, 5.6610441994727472, + 10.87954812525739}, + HypoTestInverterParams{"Hybrid_ProfileLR_Obs10_30", kHybrid, kProfileLR, 10, + 30, 7.1079825543826782, 12.553916197221213})), + backendParamName); + +/////////////////////////////////////////////////////////////////////////////// +// +// HYPOTESTINVERTER UPPER LIMIT - SIGNAL + BACKGROUND + EFFICIENCY MODEL +// +// Test the validity of the upper limit computed by the HypoTestInverter +// on a complex model distribution with signal, background and efficiency. +// Reference values and test values are both computed with the +// HypoTestInverter. As such, this test can only confirm if the +// HypoTestInverter has the same behaviour across different computing platforms +// or RooStats revisions. +// +// ModelConfig (explicit) : Poisson Signal + Background + Efficiency +// built in stressRooStats_models.h +// +// Original tests 44-48 (TestHypoTestInverter2). +// +/////////////////////////////////////////////////////////////////////////////// + +struct HypoTestInverterUpperLimitParams { + std::string name; + ECalculatorType calculatorType; + ETestStatType testStatType; + int obsValueX; // observed value "x" when measuring sig * eff + bkg + double confidenceLevel; + double refUpperLimit; + double refExpUpperLimit; + double refExpUpperLimitMinus2; + double refExpUpperLimitMinus1; + double refExpUpperLimitPlus1; + double refExpUpperLimitPlus2; +}; + +class HypoTestInverterUpperLimit : public StressRooStatsParamTest {}; + +TEST_P(HypoTestInverterUpperLimit, RegressionUpperLimit) { - bool doWrite = false; - int verbose = 0; - bool allTests = false; - bool oneTest = false; - int testNumber = 0; - bool dryRun = false; - bool doDump = false; - bool doTreeStore = false; - string refFileName = "stressRooStats_ref.root"; - - // in interpreted mode, the minimizer is hardcoded to Minuit 1 - ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit"); - - return stressRooStats(refFileName.c_str(), doWrite, verbose, allTests, oneTest, testNumber, dryRun, doDump, - doTreeStore); + // larger value test tolerance especially when using toys (difference of + // <~ 0.1 observed between using Minuit or Minuit2) + const double tol = (param().calculatorType == kAsymptotic) ? 0.02 : 0.1; + + // Create workspace and model + RooWorkspace ws{"w"}; + buildPoissonEfficiencyModel(ws); + auto sbModel = dynamic_cast(ws.obj("S+B")); + auto bModel = dynamic_cast(ws.obj("B")); + ASSERT_NE(sbModel, nullptr); + ASSERT_NE(bModel, nullptr); + + // add observed values to data set + ws.var("x")->setVal(param().obsValueX); + ws.data("data")->add(*sbModel->GetObservables()); + + // set snapshots + sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); + ws.var("sig")->setVal(0); + bModel->SetSnapshot(*bModel->GetParametersOfInterest()); + + // calculate upper limit with HypoTestInverter (the test statistic is + // declared first because the calculator will hold a raw pointer to it) + std::unique_ptr testStat{buildTestStatistic(param().testStatType, *sbModel, *bModel)}; + std::unique_ptr calc{ + buildHypoTestCalculator(param().calculatorType, *ws.data("data"), *sbModel, *bModel, 100, 100)}; + HypoTestInverter hti{*calc, nullptr, 1.0 - param().confidenceLevel}; + hti.SetTestStatistic(*testStat); + + int nscanPoints = 10; + if (param().calculatorType == kAsymptotic) { + static_cast(*calc).SetOneSided(true); + nscanPoints = 40; + } + + hti.SetFixedScan(nscanPoints, ws.var("sig")->getMin(), ws.var("sig")->getMax()); // significant speedup + + // needed because we have no extended pdf and the ToyMC Sampler evaluation + // returns an error + auto tmcs = static_cast(hti.GetHypoTestCalculator()->GetTestStatSampler()); + tmcs->SetNEventsPerToy(1); + tmcs->SetUseMultiGen(true); // make ToyMCSampler faster + + // calculate interval and extract observed upper limit and expected upper + // limit (+- 1, 2 sigma) + std::unique_ptr interval{hti.GetInterval()}; + EXPECT_NEAR(interval->UpperLimit(), param().refUpperLimit, tol); + EXPECT_NEAR(interval->GetExpectedUpperLimit(0), param().refExpUpperLimit, tol); + EXPECT_NEAR(interval->GetExpectedUpperLimit(-2), param().refExpUpperLimitMinus2, tol); + EXPECT_NEAR(interval->GetExpectedUpperLimit(-1), param().refExpUpperLimitMinus1, tol); + EXPECT_NEAR(interval->GetExpectedUpperLimit(1), param().refExpUpperLimitPlus1, tol); + EXPECT_NEAR(interval->GetExpectedUpperLimit(2), param().refExpUpperLimitPlus2, tol); } -#endif +INSTANTIATE_TEST_SUITE_P( + RooStats, HypoTestInverterUpperLimit, + ::testing::Combine(::testing::Values(ROOFIT_EVAL_BACKENDS), + ::testing::Values(HypoTestInverterUpperLimitParams{"Asymptotic_ProfileLROneSided_Obs10_CL0p95", + kAsymptotic, kProfileLROneSided, 10, 0.95, + 24.415175188860971, 11.71352859205691, 0.0, + 4.1446917240185863, 21.580917321964726, + 34.586630899315331}, + HypoTestInverterUpperLimitParams{"Asymptotic_ProfileLROneSided_Obs20_CL1Sigma", + kAsymptotic, kProfileLROneSided, 20, + kCL1Sigma, 35.455729668223988, + 3.8701997338910363, 0.0, 0.0, + 13.115199373487766, 24.476388830684247}, + HypoTestInverterUpperLimitParams{"Frequentist_RatioLR_Obs10_CL0p95", + kFrequentist, kRatioLR, 10, 0.95, + 27.777734839933814, 13.425920689335342, + 5.2777802875157622, 6.4814835475786641, + 20.000008863407345, 30.555555555555557}, + HypoTestInverterUpperLimitParams{"Frequentist_ProfileLROneSided_Obs10_CL0p95", + kFrequentist, kProfileLROneSided, 10, 0.95, + 23.611111111111089, 14.444453307851788, + 5.2777802875157622, 6.4814835475786641, + 20.251355120121389, 29.629641114958538}, + HypoTestInverterUpperLimitParams{"Hybrid_SimpleLR_Obs10_CL0p95", kHybrid, + kSimpleLR, 10, 0.95, 25.0, + 11.111095211258252, 5.2777802875157622, + 7.6388877462703384, 20.201996699382633, + 26.376273186118574})), + backendParamName); diff --git a/roofit/roostats/test/stressRooStats_ref.root b/roofit/roostats/test/stressRooStats_ref.root deleted file mode 100644 index ad3cf3f00c295..0000000000000 Binary files a/roofit/roostats/test/stressRooStats_ref.root and /dev/null differ diff --git a/roofit/roostats/test/stressRooStats_tests.h b/roofit/roostats/test/stressRooStats_tests.h deleted file mode 100644 index 9a958201211f9..0000000000000 --- a/roofit/roostats/test/stressRooStats_tests.h +++ /dev/null @@ -1,1792 +0,0 @@ -#include "stressRooStats_models.h" // Global functions that build complex RooStats models - -// RooStats headers -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -// ROOT headers -#include -#include - -// RooFit headers -#include -#include -#include -#include -#include -#include - -// STL headers -#include - -using namespace ROOT::Math; -using namespace RooFit; -using namespace RooStats; - -// testStatType = 0 Simple Likelihood Ratio (the LEP TestStat) -// = 1 Ratio of Profiled Likelihood Ratios (the Tevatron TestStat) -// = 2 Profile Likelihood Ratio (the LHC TestStat) -// = 3 Profile Likelihood One Sided (pll = 0 if mu < mu_hat) -// = 4 Profile Likelihood Signed (pll = -pll if mu < mu_hat) -// = 5 Max Likelihood Estimate as test statistic -// = 6 Number of Observed Events as test statistic -enum ECalculatorType { kAsymptotic = 0, kFrequentist = 1, kHybrid = 2 }; -enum ETestStatType { - kSimpleLR = 0, - kRatioLR = 1, - kProfileLR = 2, - kProfileLROneSided = 3, - kProfileLROneSidedDiscovery = 4, - kProfileLRSigned = 5, - kMLE = 6, - kNObs = 7 -}; -static const char *const kECalculatorTypeString[] = {"Asymptotic", "Frequentist", "Hybrid"}; -static const char *const kETestStatTypeString[] = {"Simple-Likelihood-Ratio", - "Ratio-Of-Profiled-Likelihoods", - "Profile-Likelihood-Ratio", - "Profile-Likelihood-One-Sided", - "Profile-Likelihood-One-Sided-Discovery", - "Profile-Likelihood-Signed", - "Max-Likelihood-Estimate", - "Number-Of-Observed-Events"}; -// static const char * const kETestStatTypeString[] = { "Simple Likelihood Ratio", "Ratio Of Profiled Likelihoods", -// "Profile Likelihood Ratio", "Profile Likelihood One-Sided", "Profile Likelihood One-Sided Discovery", -// "Profile Likelihood Signed", "Max Likelihood Estimate", "Number Of Observed Events" }; -static HypoTestCalculatorGeneric *buildHypoTestCalculator(const ECalculatorType calculatorType, RooAbsData &data, - const ModelConfig &nullModel, const ModelConfig &altModel, - const UInt_t toysNull, const UInt_t toysAlt); -static TestStatistic * -buildTestStatistic(const ETestStatType testStatType, const ModelConfig &sbModel, const ModelConfig &bModel); - -//_____________________________________________________________________________ -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -// -// PART ONE: -// PROFILE LIKELIHOOD CALCULATOR UNIT TESTS -// - -/////////////////////////////////////////////////////////////////////////////// -// -// PROFILE LIKELIHOOD CALCULATOR - LIKELIHOOD INTERVAL - GAUSSIAN DISTRIBUTION -// -// Test the likelihood interval computed by the profile likelihood calculator -// on a Gaussian distribution. Reference interval limits are computed via -// analytic methods: solve equation 2*(ln(LL(xMax))-ln(LL(x)) = q, where q = -// normal_quantile_c(testSize/2, 1). In the case of a Gaussian distribution, the -// interval limits are equal to: mean +- normal_quantile_c(testSize/2, sigma/sqrt(N)). -// -// ModelConfig (implicit) : -// Observable -> x -// Parameter of Interest -> mean -// Nuisance parameter (Constant !) -> sigma -// -// Input Parameters: -// confidenceLevel -> Confidence Level of the interval we are calculating -// -// 03/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestProfileLikelihoodCalculator1 : public RooUnitTest { -private: - double fConfidenceLevel; - -public: - TestProfileLikelihoodCalculator1(TFile *refFile, bool writeRef, Int_t verbose, double confidenceLevel = 0.95) - : RooUnitTest("ProfileLikelihoodCalculator Interval - Gaussian Model", refFile, writeRef, verbose), - fConfidenceLevel(confidenceLevel){}; - - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fConfidenceLevel <= 0.0 || fConfidenceLevel >= 1.0) { - Warning("isTestAvailable", "Confidence level must be in the range (0,1). Skipping test..."); - return false; - } - return true; - } - - bool testCode() override - { - - const Int_t N = 10; // number of observations - // the compared values / objects must have the same name in write / compare modes - const TString lowerLimitString = TString::Format("tplc2_lower_limit_mean_%lf", fConfidenceLevel); - const TString upperLimitString = TString::Format("tplc2_upper_limit_mean_%lf", fConfidenceLevel); - - // TODO: see why it fails for a small number of observations - // Create Gaussian model, generate data set and define - RooWorkspace ws{"w"}; - ws.factory("Gaussian::gauss(x[-5,5], mean[0,-5,5], sigma[1])"); - std::unique_ptr data{ws.pdf("gauss")->generate(*ws.var("x"), N)}; - - if (_write == true) { - - // Calculate likelihood interval from data via analytic methods - double estMean = data->mean(*ws.var("x")); - double intervalHalfWidth = - normal_quantile_c((1.0 - fConfidenceLevel) / 2.0, ws.var("sigma")->getValV() / sqrt((double)N)); - double lowerLimit = estMean - intervalHalfWidth; - double upperLimit = estMean + intervalHalfWidth; - - // Compare the limits obtained via ProfileLikelihoodCalculator with analytically estimated values - regValue(lowerLimit, lowerLimitString); - regValue(upperLimit, upperLimitString); - - } else { - - // Calculate likelihood interval using the ProfileLikelihoodCalculator - auto plc = std::make_unique(*data, *ws.pdf("gauss"), *ws.var("mean")); - plc->SetConfidenceLevel(fConfidenceLevel); - std::unique_ptr interval{plc->GetInterval()}; - - // Register analytically computed limits in the reference file - regValue(interval->LowerLimit(*ws.var("mean")), lowerLimitString); - regValue(interval->UpperLimit(*ws.var("mean")), upperLimitString); - - plc.reset(); - } - - return true; - } -}; - -/////////////////////////////////////////////////////////////////////////////// -// -// PROFILE LIKELIHOOD CALCULATOR - LIKELIHOOD INTERVAL - POISSON DISTRIBUTION -// -// Test the 68% likelihood interval computed by the profile likelihood calculator -// on a Poisson distribution, from only one observed value. Reference values are -// computed via analytic methods: solve equation 2*[ln(LL(xMax)) - ln(LL(x))] = 1. -// -// ModelConfig (implicit) : -// Observable -> x -// Parameter of Interest -> mean -// -// Input Parameters: -// obsValue -> observed value in experiment -// -// 03/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestProfileLikelihoodCalculator2 : public RooUnitTest { -private: - Int_t fObsValue; - -public: - TestProfileLikelihoodCalculator2(TFile *refFile, bool writeRef, Int_t verbose, Int_t obsValue = 5) - : RooUnitTest("ProfileLikelihoodCalculator Interval - Poisson Simple Model", refFile, writeRef, verbose), - fObsValue(obsValue){}; - - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fObsValue < 0 || fObsValue > 1000) { - Warning("isTestAvailable", "Observed value must be in the range [0,1000]. Skipping test..."); - return false; - } - return true; - } - - bool testCode() override - { - - // the compared values / objects must have the same name in write / compare modes - const TString lowerLimitString = TString::Format("tplc2_lower_limit_mean_%d", fObsValue); - const TString upperLimitString = TString::Format("tplc2_upper_limit_mean_%d", fObsValue); - - // write reference values - if (_write == true) { - - // Solutions of equation 2*[ln(LL(xMax)) - ln(LL(x))] = 1, where xMax is the point of maximum likelihood - // For the special case of the Poisson distribution with N = 1, xMax = obsValue - TString llRatioExpression = - TString::Format("2*(x-%d*log(x)-%d+%d*log(%d))", fObsValue, fObsValue, fObsValue, fObsValue); - // Special case fObsValue = 0 because log(0) not computable, the limit of n * log(n), n->0 must be taken - if (fObsValue == 0) - llRatioExpression = TString::Format("2*x"); - - auto llRatio = std::make_unique("llRatio", llRatioExpression, 1e-100, fObsValue); // lowerLimit < obsValue - double lowerLimit = llRatio->GetX(1); - llRatio->SetRange(fObsValue, 1000); // upperLimit > obsValue - double upperLimit = llRatio->GetX(1); - - // Compare the limits obtained via ProfileLikelihoodCalculator with the likelihood ratio analytic computations - regValue(lowerLimit, lowerLimitString); - regValue(upperLimit, upperLimitString); - - // compare with reference values - } else { - - // Set a 68% confidence level for the interval - const double confidenceLevel = 2 * normal_cdf(1) - 1.0; - - // Create Poisson model and dataset - RooWorkspace ws{"w"}; - ws.factory(TString::Format("Poisson::poiss(x[%d,0,1000], mean[0,1000])", fObsValue).Data()); - RooDataSet data{"data", "data", *ws.var("x")}; - data.add(*ws.var("x")); - - // Calculate likelihood interval using the ProfileLikelihoodCalculator - auto plc = std::make_unique(data, *ws.pdf("poiss"), *ws.var("mean")); - plc->SetConfidenceLevel(confidenceLevel); - std::unique_ptr interval{plc->GetInterval()}; - - // Register externally computed limits in the reference file - regValue(interval->LowerLimit(*ws.var("mean")), lowerLimitString); - regValue(interval->UpperLimit(*ws.var("mean")), upperLimitString); - - plc.reset(); - } - - return true; - } -}; - -/////////////////////////////////////////////////////////////////////////////// -// -// PROFILE LIKELIHOOD CALCULATOR - LIKELIHOOD INTERVAL - POISSON PRODUCT MODEL -// -// Test the 68% likelihood interval computed by the ProfileLikelihoodCalculator -// on a complex model. Reference values and test values are both computed with -// the ProfileLikelihoodCalculator. As such, this test can only confirm if the -// ProfileLikelihoodCalculator has the same behaviour across different computer -// platforms or RooStats revisions. -// -// ModelConfig (explicit) : Poisson Product Model -// built in stressRooStats_models.cxx -// -// Input Parameters: -// obsValueX -> observed value "x" when measuring sig + bkg1 -// obsValueY -> observed value "y" when measuring 2*sig*1.2^beta + bkg2 -// confidenceLevel -> Confidence Level of the interval we are calculating -// -// 04/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestProfileLikelihoodCalculator3 : public RooUnitTest { -private: - Int_t fObsValueX; - Int_t fObsValueY; - double fConfidenceLevel; - -public: - TestProfileLikelihoodCalculator3(TFile *refFile, bool writeRef, Int_t verbose, Int_t obsValueX = 15, - Int_t obsValueY = 30, double confidenceLevel = 2 * normal_cdf(1) - 1) - : RooUnitTest("ProfileLikelihoodCalculator Interval - Poisson Product Model", refFile, writeRef, verbose), - fObsValueX(obsValueX), - fObsValueY(obsValueY), - fConfidenceLevel(confidenceLevel){}; - - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fObsValueX < 0 || fObsValueX > 30) { - Warning("isTestAvailable", "Observed value X=s+b must be in the range [0,30]. Skipping test..."); - return false; - } - if (fObsValueY < 0 || fObsValueY > 80) { - Warning("isTestAvailable", "Observed value Y=2*s*1.2^beta+b must be in the range [0,80]. Skipping test..."); - return false; - } - if (fConfidenceLevel <= 0.0 || fConfidenceLevel >= 1.0) { - Warning("isTestAvailable", "Confidence level must be in the range (0,1). Skipping test..."); - return false; - } - return true; - } - - bool testCode() override - { - - // Create workspace and model - RooWorkspace ws{"w"}; - buildPoissonProductModel(&ws); - auto model = dynamic_cast(ws.obj("S+B")); - - // add observed values to data set - ws.var("x")->setVal(fObsValueX); - ws.var("y")->setVal(fObsValueY); - ws.data("data")->add(*model->GetObservables()); - - std::unique_ptr initialVariables{model->GetPdf()->getVariables()}; - ws.saveSnapshot("initialVariables", *initialVariables); - - // build likelihood interval with ProfileLikelihoodCalculator - ProfileLikelihoodCalculator plc{*ws.data("data"), *model}; - plc.SetConfidenceLevel(fConfidenceLevel); - std::unique_ptr interval{plc.GetInterval()}; - - regValue(interval->LowerLimit(*ws.var("sig")), - TString::Format("tplc3_lower_limit_sig_%d_%d_%lf", fObsValueX, fObsValueY, fConfidenceLevel)); - regValue(interval->UpperLimit(*ws.var("sig")), - TString::Format("tplc3_upper_limit_sig_%d_%d_%lf", fObsValueX, fObsValueY, fConfidenceLevel)); - - if (_verb > 1) { - ws.loadSnapshot("initialVariables"); - ws.writeToFile(TString::Format("stressRooStats_PoissonProductModel_%d_%d.root", fObsValueX, fObsValueY)); - } - - return true; - } -}; - -/////////////////////////////////////////////////////////////////////////////// -// -// PROFILE LIKELIHOOD CALCULATOR HYPOTHESIS TEST - ON / OFF MODEL -// -// Perform a hypothesis test using the ProfileLikelihoodCalculator on the -// on/off model. Reference values and test values are both computed with the -// ProfileLikelihoodCalculator. As such, this test can only confirm if the -// ProfileLikelihoodCalculator has the same behaviour across different -// computing platforms or RooStats revisions. -// -// ModelConfig (explicit) : Poisson On / Off Model -// built in stressRooStats_models.cxx -// -// For a detailed description of the on/off model, see the paper: "Evaluation -// of three methods for calculating statistical significance when incorporating -// a systematic uncertainty into a test of the background-only hypothesis for -// a Poisson process" by Robert D. Cousins, James T. Linnemann, Jordan Tucker -// -// 04/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestProfileLikelihoodCalculator4 : public RooUnitTest { -public: - TestProfileLikelihoodCalculator4(TFile *refFile, bool writeRef, Int_t verbose) - : RooUnitTest("ProfileLikelihoodCalculator Hypothesis Test", refFile, writeRef, verbose){}; - - // Override test value tolerance - // A larger tolerance is needed since the values in the Cousins paper are given with 1e-2 precision - double vtol() override { return 1e-2; } - - bool testCode() override - { - - // For testing purposes, we consider four special cases for which the values are known from - // the Cousins et al. paper mentioned above. The inputs for each of these cases are (using - // the notations from the paper): n_on, n_off and Z_PL. We provide a certain fixed input set - // for each case. - const Int_t numberTestSets = 3; - const Int_t numberOnEvents[numberTestSets] = {4, 50, 67}; - const Int_t numberOffEvents[numberTestSets] = {5, 55, 15}; - const double tau[numberTestSets] = {5.0, 2.0, 0.5}; - const double significance[numberTestSets] = {1.95, 3.02, 3.04}; - - for (Int_t i = 0; i < numberTestSets; ++i) { - - TString stringSignificance = - TString::Format("tplc4_significance_%d_%d_%lf", numberOnEvents[i], numberOffEvents[i], tau[i]); - - if (_write == true) { - - // register reference values from Cousins et al. paper - regValue(significance[i], stringSignificance); - - } else { - - // build workspace and model - RooWorkspace ws{"w"}; - buildOnOffModel(ws); - auto sbModel = dynamic_cast(ws.obj("S+B")); - auto bModel = dynamic_cast(ws.obj("B")); - - // add observable values to data set - ws.var("n_on")->setVal(numberOnEvents[i]); - ws.var("n_off")->setVal(numberOffEvents[i]); - ws.var("tau")->setVal(tau[i]); - ws.var("tau")->setConstant(); - ws.data("data")->add(*sbModel->GetObservables()); - - // set snapshots - ws.var("sig")->setVal(numberOnEvents[i] - numberOffEvents[i] / tau[i]); - sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); - ws.var("sig")->setVal(0); - bModel->SetSnapshot(*bModel->GetParametersOfInterest()); - - // has as initial value a non-zero value for sig (i.e start with the S+B value) - sbModel->LoadSnapshot(); - - // get significance using the ProfileLikelihoodCalculator - auto plc = std::make_unique(*ws.data("data"), *sbModel); - plc->SetNullParameters(*bModel->GetSnapshot()); - // plc->SetAlternateParameters(*sbModel->GetSnapshot()); // not needed for PLC - - regValue(plc->GetHypoTest()->Significance(), stringSignificance); - } - } - - return true; - } -}; - -// -// END OF PART ONE -// -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -//_____________________________________________________________________________ - -//_____________________________________________________________________________ -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -// -// PART TWO: -// BAYESIAN CALCULATOR UNIT TESTS -// - -/////////////////////////////////////////////////////////////////////////////// -// -// BAYESIAN CENTRAL INTERVAL - SIMPLE MODEL -// -// Test the Bayesian central interval computed by the BayesianCalculator on a -// Poisson distribution, using different priors. The parameter of interest is -// the mean of the Poisson distribution, and there are no nuisance parameters. -// The priors used are: -// 1. constant / uniform -// 2. inverse of the mean -// 3. square root of the inverse of the mean -// 4. gamma distribution -// The posterior distribution is easily obtained analytically for these cases. -// Therefore, the reference interval limits will be computed analytically. -// -// ModelConfig (implicit) : -// Observable -> x -// Parameter of Interest -> mean -// -// Input Parameters: -// obsValue -> observed value in experiment -// confidenceLevel -> Confidence Level of the interval we are calculating -// -// 04/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestBayesianCalculator1 : public RooUnitTest { -private: - Int_t fObsValue; - double fConfidenceLevel; - static double priorInv(double mean) { return 1.0 / mean; } - static double priorInvSqrt(double mean) { return 1.0 / sqrt(mean); } - -public: - TestBayesianCalculator1(TFile *refFile, bool writeRef, Int_t verbose, Int_t obsValue = 3, - double confidenceLevel = 2 * normal_cdf(1) - 1) - : RooUnitTest("BayesianCalculator Central Interval - Poisson Simple Model", refFile, writeRef, verbose), - fObsValue(obsValue), - fConfidenceLevel(confidenceLevel){}; - - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fObsValue < 0 || fObsValue > 100) { - Warning("isTestAvailable", "Observed value must be in the range [0,100]. Skipping test..."); - return false; - } - if (fConfidenceLevel <= 0.0 || fConfidenceLevel >= 1.0) { - Warning("isTestAvailable", "Confidence level must be in the range (0,1). Skipping test..."); - return false; - } - return true; - } - - bool testCode() override - { - - // Set the confidence level for a 68.3% CL central interval - const double gammaShape = 2; // shape of the gamma distribution prior (gamma = alpha) - const double gammaRate = 1; // rate = 1/scale of the gamma distribution prior (beta = 1/theta) - const Int_t numberScans = 10000; // tested to be sufficient for the scan of the Bayesian posterior - - // names of tested variables must be the same in write / comparison modes - const TString lowerLimitString = TString::Format("tbc1_lower_limit_unif_%d_%lf", fObsValue, fConfidenceLevel); - const TString upperLimitString = TString::Format("tbc1_upper_limit_unif_%d_%lf", fObsValue, fConfidenceLevel); - const TString lowerLimitInvString = TString::Format("tbc1_lower_limit_inv_%d_%lf", fObsValue, fConfidenceLevel); - const TString upperLimitInvString = TString::Format("tbc1_upper_limit_inv_%d_%lf", fObsValue, fConfidenceLevel); - const TString lowerLimitInvSqrtString = - TString::Format("tbc1_lower_limit_inv_sqrt_%d_%lf", fObsValue, fConfidenceLevel); - const TString upperLimitInvSqrtString = - TString::Format("tbc1_upper_limit_inv_sqrt_%d_%lf", fObsValue, fConfidenceLevel); - const TString lowerLimitGammaString = - TString::Format("tbc1_lower_limit_gamma_%d_%lf", fObsValue, fConfidenceLevel); - const TString upperLimitGammaString = - TString::Format("tbc1_upper_limit_gamma_%d_%lf", fObsValue, fConfidenceLevel); - - if (_write == true) { - - double lowerLimit = gamma_quantile((1.0 - fConfidenceLevel) / 2, fObsValue + 1, 1); // integrate to 16% - double upperLimit = gamma_quantile_c((1.0 - fConfidenceLevel) / 2, fObsValue + 1, 1); // integrate to 84% - double lowerLimitInv = gamma_quantile((1.0 - fConfidenceLevel) / 2, fObsValue, 1); - double upperLimitInv = gamma_quantile_c((1.0 - fConfidenceLevel) / 2, fObsValue, 1); - double lowerLimitInvSqrt = gamma_quantile((1.0 - fConfidenceLevel) / 2, fObsValue + 0.5, 1); - double upperLimitInvSqrt = gamma_quantile_c((1.0 - fConfidenceLevel) / 2, fObsValue + 0.5, 1); - double lowerLimitGamma = - gamma_quantile((1.0 - fConfidenceLevel) / 2, fObsValue + gammaShape, 1.0 / (1 + gammaRate)); - double upperLimitGamma = - gamma_quantile_c((1.0 - fConfidenceLevel) / 2, fObsValue + gammaShape, 1.0 / (1 + gammaRate)); - - // Compare the limits obtained via BayesianCalculator with quantile values - regValue(lowerLimit, lowerLimitString); - regValue(upperLimit, upperLimitString); - regValue(lowerLimitInv, lowerLimitInvString); - regValue(upperLimitInv, upperLimitInvString); - regValue(lowerLimitInvSqrt, lowerLimitInvSqrtString); - regValue(upperLimitInvSqrt, upperLimitInvSqrtString); - regValue(lowerLimitGamma, lowerLimitGammaString); - regValue(upperLimitGamma, upperLimitGammaString); - - } else { - - // Create Poisson model - RooWorkspace ws{"w"}; - ws.factory("Poisson::poiss(x[0,100], mean[1e-6,100])"); - // TODO: see why it does not work so well for boundary observed values {0, 100} - - // create prior pdfs - ws.factory("Uniform::prior(mean)"); - ws.import(RooCFunction1PdfBinding("priorInv", "priorInv", &priorInv, *ws.var("mean"))); - ws.import( - RooCFunction1PdfBinding("priorInvSqrt", "priorInvSqrt", priorInvSqrt, *ws.var("mean"))); - ws.factory(TString::Format("Gamma::priorGamma(mean, %lf, %lf, 0)", gammaShape, gammaRate).Data()); - - // build argument sets and data set - ws.defineSet("obs", "x"); - ws.defineSet("poi", "mean"); - ws.var("x")->setVal(fObsValue); - RooDataSet data{"data", "data", *ws.set("obs")}; - data.add(*ws.set("obs")); - - // NOTE: RooIntegrator1D is too slow and gives poor results -#ifdef ROOFITMORE - RooAbsReal::defaultIntegratorConfig()->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D"); -#endif - - std::unique_ptr bc; - std::unique_ptr interval; - - // Uniform prior on mean - bc = std::make_unique(data, *ws.pdf("poiss"), *ws.set("poi"), *ws.pdf("prior"), nullptr); - bc->SetConfidenceLevel(fConfidenceLevel); - bc->SetScanOfPosterior(numberScans); - interval = std::unique_ptr{bc->GetInterval()}; - regValue(interval->LowerLimit(), lowerLimitString); - regValue(interval->UpperLimit(), upperLimitString); - - // Inverse of mean prior - bc = - std::make_unique(data, *ws.pdf("poiss"), *ws.set("poi"), *ws.pdf("priorInv"), nullptr); - bc->SetConfidenceLevel(fConfidenceLevel); - bc->SetScanOfPosterior(numberScans); - interval = std::unique_ptr{bc->GetInterval()}; - regValue(interval->LowerLimit(), lowerLimitInvString); - regValue(interval->UpperLimit(), upperLimitInvString); - - // Square root of inverse of mean prior - bc = std::make_unique(data, *ws.pdf("poiss"), *ws.set("poi"), *ws.pdf("priorInvSqrt"), - nullptr); - bc->SetConfidenceLevel(fConfidenceLevel); - bc->SetScanOfPosterior(numberScans); - interval = std::unique_ptr{bc->GetInterval()}; - regValue(interval->LowerLimit(), lowerLimitInvSqrtString); - regValue(interval->UpperLimit(), upperLimitInvSqrtString); - - // Gamma distribution prior - bc = std::make_unique(data, *ws.pdf("poiss"), *ws.set("poi"), *ws.pdf("priorGamma"), - nullptr); - bc->SetConfidenceLevel(fConfidenceLevel); - bc->SetScanOfPosterior(numberScans); - interval = std::unique_ptr{bc->GetInterval()}; - regValue(interval->LowerLimit(), lowerLimitGammaString); - regValue(interval->UpperLimit(), upperLimitGammaString); - } - - return true; - } -}; - -/////////////////////////////////////////////////////////////////////////////// -// -// BAYESIAN SHORTEST INTERVAL - SIMPLE POISSON MODEL -// -// Test the Bayesian shortest interval computed by the BayesianCalculator on a -// Poisson distribution, using different priors. The parameter of interest is -// the mean of the Poisson distribution, and there are no nuisance parameters. -// The priors used are: -// 1. constant / uniform -// 2. inverse of the mean -// The reference interval limits are taken from the paper: "Why isn't every -// physicist a Bayesian?" by Robert D. Cousins. -// -// ModelConfig (implicit) : -// Observable -> x -// Parameter of Interest -> mean -// -// 04/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestBayesianCalculator2 : public RooUnitTest { -public: - TestBayesianCalculator2(TFile *refFile, bool writeRef, Int_t verbose) - : RooUnitTest("BayesianCalculator Shortest Interval - Poisson Simple Model", refFile, writeRef, verbose){}; - - // the references values in the paper have a precision of only two decimal points - // in such a situation, it is natural that we increase the value tolerance - double vtol() override { return 1e-2; } - - bool testCode() override - { - - // Put the confidence level so that we obtain a 68% confidence interval - const double confidenceLevel = 2 * normal_cdf(1) - 1; - const Int_t obsValue = 3; // observed experiment value - const Int_t numberScans = 100000; // sufficient number of scans - - // names of tested variables must be the same in write / comparison modes - const TString lowerLimitString = "tbc2_lower_limit_unif"; - const TString upperLimitString = "tbc2_upper_limit_unif"; - const TString lowerLimitInvString = "tbc2_lower_limit_inv"; - const TString upperLimitInvString = "tbc2_upper_limit_inv"; - - if (_write == true) { - - // Compare the limits obtained via BayesianCalculator with given reference values - regValue(1.55, lowerLimitString); - regValue(5.15, upperLimitString); - regValue(0.86, lowerLimitInvString); - regValue(3.85, upperLimitInvString); - - } else { - - // Create Poisson model - RooWorkspace ws{"w"}; - ws.factory("Poisson::poiss(x[0,100], mean[1e-6,100])"); - ws.factory("Uniform::prior(mean)"); - ws.factory("EXPR::priorInv('1/mean', mean)"); - - // build argument sets and data set - ws.defineSet("poi", "mean"); - ws.defineSet("obs", "x"); - ws.var("x")->setVal(obsValue); - RooDataSet data{"data", "data", *ws.set("obs")}; - data.add(*ws.set("obs")); - - // NOTE: RooIntegrator1D is too slow and gives poor results -#ifdef ROOFITMORE - RooAbsReal::defaultIntegratorConfig()->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D"); -#endif - // Uniform prior on mean - auto bc = std::make_unique(data, *ws.pdf("poiss"), *ws.set("poi"), *ws.pdf("prior"), nullptr); - bc->SetConfidenceLevel(confidenceLevel); - bc->SetShortestInterval(); - bc->SetScanOfPosterior(numberScans); - std::unique_ptr interval{bc->GetInterval()}; - regValue(interval->LowerLimit(), lowerLimitString); - regValue(interval->UpperLimit(), upperLimitString); - - // Inverse of mean prior - bc = std::make_unique(data, *ws.pdf("poiss"), *ws.set("poi"), *ws.pdf("priorInv"), nullptr); - bc->SetConfidenceLevel(confidenceLevel); - bc->SetShortestInterval(); - bc->SetScanOfPosterior(numberScans); - interval = std::unique_ptr{bc->GetInterval()}; - regValue(interval->LowerLimit(), lowerLimitInvString); - regValue(interval->UpperLimit(), upperLimitInvString); - } - - return true; - } -}; - -/////////////////////////////////////////////////////////////////////////////// -// -// BAYESIAN CENTRAL INTERVAL - POISSON PRODUCT MODEL -// -// Test the validity of the central interval computed by the BayesianCalculator -// on a complex Poisson model distribution. Reference values and test values -// are both computed with the BayesianCalculator. As such, this test can only -// confirm if the BayesianCalculator has the same behaviour across different -// computing platforms or RooStats revisions. A uniform prior PDF is used for the -// parameter of interest ("sig"). -// -// ModelConfig (explicit) : Poisson Product Model -// built in stressRooStats_models.cxx -// -// Input Parameters: -// obsValueX -> observed value "x" when measuring sig + bkg1 -// obsValueY -> observed value "y" when measuring 2*sig*1.2^beta + bkg2 -// confidenceLevel -> Confidence Level of the interval we are calculating -// -// 04/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestBayesianCalculator3 : public RooUnitTest { -private: - Int_t fObsValueX; - Int_t fObsValueY; - double fConfidenceLevel; - -public: - TestBayesianCalculator3(TFile *refFile, bool writeRef, Int_t verbose, Int_t obsValueX = 15, Int_t obsValueY = 30, - double confidenceLevel = 2 * normal_cdf(1) - 1) - : RooUnitTest("BayesianCalculator Central Interval - Poisson Product Model", refFile, writeRef, verbose), - fObsValueX(obsValueX), - fObsValueY(obsValueY), - fConfidenceLevel(confidenceLevel){}; - - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fObsValueX < 0 || fObsValueX > 30) { - Warning("isTestAvailable", "Observed value X=s+b must be in the range [0,30]. Skipping test..."); - return false; - } - if (fObsValueY < 0 || fObsValueY > 80) { - Warning("isTestAvailable", "Observed value Y=2*s*1.2^beta+b must be in the range [0,80]. Skipping test..."); - return false; - } - if (fConfidenceLevel <= 0.0 || fConfidenceLevel >= 1.0) { - Warning("isTestAvailable", "Confidence level must be in the range (0,1). Skipping test..."); - return false; - } - return true; - } - - bool testCode() override - { - - const Int_t numberScans = 10; // sufficient number of scans - - // Create workspace and model - auto w = std::make_unique("w"); - buildPoissonProductModel(w.get()); - ModelConfig *model = static_cast(w->obj("S+B")); - - // add observed values to data set - w->var("x")->setVal(fObsValueX); - w->var("y")->setVal(fObsValueY); - w->data("data")->add(*model->GetObservables()); - - std::unique_ptr initialVariables{model->GetPdf()->getVariables()}; - w->saveSnapshot("initialVariables", *initialVariables); - - // NOTE: Roo1DIntegrator is too slow and gives poor results -#ifdef ROOFITMORE - RooAbsReal::defaultIntegratorConfig()->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D"); -#endif - - // Create BayesianCalculator and - auto bc = std::make_unique(*w->data("data"), *model); - bc->SetConfidenceLevel(fConfidenceLevel); - bc->SetScanOfPosterior(numberScans); - - // Obtain confidence interval by scanning the posterior function in the given number of points - std::unique_ptr interval{bc->GetInterval()}; - regValue(interval->LowerLimit(), - TString::Format("tbc3_lower_limit_sig_%d_%d_%lf", fObsValueX, fObsValueY, fConfidenceLevel)); - regValue(interval->UpperLimit(), - TString::Format("tbc3_upper_limit_sig_%d_%d_%lf", fObsValueX, fObsValueY, fConfidenceLevel)); - - return true; - } -}; - -// -// END OF PART TWO -// -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -//_____________________________________________________________________________ - -//_____________________________________________________________________________ -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -// -// PART THREE: -// MARKOV CHAIN MONTE CARLO CALCULATOR UNIT TESTS -// - -/////////////////////////////////////////////////////////////////////////////// -// -// MCMC INTERVAL CALCULATOR - POISSON PRODUCT MODEL -// -// Test the validity of the confidence interval computed by the MCMCCalculator -// on a complex Poisson model distribution. Reference values and test values -// are both computed with the MCMCCalculator. As such, this test can only -// confirm if the MCMCCalculator has the same behaviour across different -// computing platforms or RooStats revisions. -// -// ModelConfig (explicit) : Poisson Product Model -// built in stressRooStats_models.cxx -// -// Input Parameters: -// obsValueX -> observed value "x" when measuring sig + bkg1 -// obsValueY -> observed value "y" when measuring 2*sig*1.2^beta + bkg2 -// confidenceLevel -> Confidence Level of the interval we are calculating -// -// 04/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestMCMCCalculator : public RooUnitTest { -private: - Int_t fObsValueX; - Int_t fObsValueY; - double fConfidenceLevel; - -public: - TestMCMCCalculator(TFile *refFile, bool writeRef, Int_t verbose, Int_t obsValueX = 15, Int_t obsValueY = 30, - double confidenceLevel = 2 * normal_cdf(1) - 1) - : RooUnitTest("MCMCCalculator Interval - Poisson Product Model", refFile, writeRef, verbose), - fObsValueX(obsValueX), - fObsValueY(obsValueY), - fConfidenceLevel(confidenceLevel){}; - - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fObsValueX < 0 || fObsValueX > 30) { - Warning("isTestAvailable", "Observed value X=s+b must be in the range [0,30]. Skipping test..."); - return false; - } - if (fObsValueY < 0 || fObsValueY > 80) { - Warning("isTestAvailable", "Observed value Y=2*s*1.2^beta+b must be in the range [0,80]. Skipping test..."); - return false; - } - if (fConfidenceLevel <= 0.0 || fConfidenceLevel >= 1.0) { - Warning("isTestAvailable", "Confidence level must be in the range (0,1). Skipping test..."); - return false; - } - return true; - } - - bool testCode() override - { - - // Create workspace and model - auto w = std::make_unique("w"); - buildPoissonProductModel(w.get()); - ModelConfig *model = static_cast(w->obj("S+B")); - - // add observed values to data set - w->var("x")->setVal(fObsValueX); - w->var("y")->setVal(fObsValueY); - w->data("data")->add(*model->GetObservables()); - - std::unique_ptr initialVariables{model->GetPdf()->getVariables()}; - w->saveSnapshot("initialVariables", *initialVariables); - - // NOTE: Roo1DIntegrator is too slow and gives poor results -#ifdef ROOFITMORE - RooAbsReal::defaultIntegratorConfig()->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D"); -#endif - - // create and configure MCMC calculator - auto sp = std::make_unique(0.1); - auto mcmcc = std::make_unique(*w->data("data"), *model); - mcmcc->SetProposalFunction(*sp); - mcmcc->SetNumIters(100000); // Metropolis-Hastings algorithm iterations - mcmcc->SetNumBurnInSteps(50); // first 50 steps to be ignored as burn-in - mcmcc->SetConfidenceLevel(fConfidenceLevel); - - // calculate the confidence interval - std::unique_ptr interval{mcmcc->GetInterval()}; - regValue(interval->LowerLimit(*w->var("sig")), - TString::Format("mcmcc_lower_limit_sig_%d_%d_%lf", fObsValueX, fObsValueY, fConfidenceLevel)); - regValue(interval->UpperLimit(*w->var("sig")), - TString::Format("mcmcc_upper_limit_sig_%d_%d_%lf", fObsValueX, fObsValueY, fConfidenceLevel)); - - return true; - } -}; - -// -// END OF PART THREE -// -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -//_____________________________________________________________________________ - -//_____________________________________________________________________________ -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -// -// PART FOUR: -// HYPOTHESIS TEST CALCULATOR UNIT TESTS -// - -///////////////////////////////////////////////////////////////////////// -// -// ZBI - ON / OFF MODEL -// -// Evaluate the functionality of the top level functions in RooStats -// called NumberCountingUtils::BinomialWithTauObsZ. This function -// computes the significance of a hypothesis test via a frequentist -// solution. This significance, called ZBi, is detailed in the article -// "Evaluation of three methods for calculating statistical significance -// when incorporating a systematic uncertainty into a test of the -// background-only hypothesis for a Poisson process" by Robert D. Cousins, -// James T. Linnemann, Jordan Tucker. The reference values are taken -// from the paper, as well as the On / Off model on which the test is -// evaluated. -// -// ModelConfig (implicit) : Poisson On / Off Model -// built in stressRooStats_models.cxx -// implicit in NumberCountingUtils::BinomialWithTauObsZ -// -// 05/2012 - Wouter Verkerke, Lorenzo Moneta, Ioan Gabriel Bucur -// -///////////////////////////////////////////////////////////////////////// - -class TestZBi : public RooUnitTest { -public: - TestZBi(TFile *refFile, bool writeRef, Int_t verbose) - : RooUnitTest("ZBi Significance - On / Off Model", refFile, writeRef, verbose){}; - - // Override test value tolerance - // A larger tolerance is needed since the values in the Cousins paper are given with 1e-2 precision - double vtol() override { return 1e-2; } - - bool testCode() override - { - - // For testing purposes, we consider four special cases for which the values are known from - // the Cousins et al. paper mentioned above. The inputs for each of these cases are (using - // the notations from the paper): n_on, n_off and Z_PL. We provide a certain fixed input set - // for each case. - const Int_t numberTestSets = 4; - const Int_t numberOnEvents[numberTestSets] = {4, 50, 67, 200}; - const Int_t numberOffEvents[numberTestSets] = {5, 55, 15, 10}; - const double tau[numberTestSets] = {5.0, 2.0, 0.5, 0.1}; - const double significance[numberTestSets] = {1.66, 2.93, 2.89, 2.2}; - - for (Int_t i = 0; i < numberTestSets; ++i) { - - TString stringSignificance = - TString::Format("tzbi_significance_%d_%d_%lf", numberOnEvents[i], numberOffEvents[i], tau[i]); - - if (_write == true) { - - // register reference values from Cousins et al. paper - regValue(significance[i], stringSignificance); - - } else { - - // call top level function - regValue(NumberCountingUtils::BinomialWithTauObsZ(numberOnEvents[i], numberOffEvents[i], tau[i]), - stringSignificance); - } - } - - return true; - } -}; - -/////////////////////////////////////////////////////////////////////////////// -// -// ASYMPTOTIC CALCULATOR VS PROFILE LIKELIHOOD CALCULATOR HYPOTHESIS TEST -// -// This test evaluates the functionality of the AsymptoticCalculator by -// comparing the significance given from a hypothesis test on the on/off model -// with the significance given by the ProfileLikelihoodCalculator. The validity -// of the PLC hypothesis test is evaluated in TestProfileLikelihoodCalculator4. -// If working properly, the two methods should yield identical results. -// -// ModelConfig (explicit) : Poisson On / Off Model -// built in stressRooStats_models.cxx -// -// Input Parameters: -// obsValueOn -> observed value "n_on" of sig + bkg -// obsValueOff -> observed value "n_off" of tau * bkg -// tau -> parameter of the model (constant with regard to integration) -// -// 05/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestHypoTestCalculator1 : public RooUnitTest { -private: - Int_t fObsValueOn; - Int_t fObsValueOff; - double fTau; - -public: - TestHypoTestCalculator1(TFile *refFile, bool writeRef, Int_t verbose, Int_t obsValueOn = 150, - Int_t obsValueOff = 100, double tau = 1.0) - : RooUnitTest("AsymptoticCalculator vs ProfileLikelihoodCalculator Significance - On / Off Model", refFile, - writeRef, verbose), - fObsValueOn(obsValueOn), - fObsValueOff(obsValueOff), - fTau(tau){}; - - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fObsValueOn < 0 || fObsValueOn > 300) { - Warning("isTestAvailable", "Observed value on_source=s+b must be in the range [0,300]. Skipping test..."); - return false; - } - if (fObsValueOff < 0 || fObsValueOff > 1100) { - Warning("isTestAvailable", "Observed value off_source=tau*b must be in the range [0,1100]. Skipping test..."); - return false; - } - if (fTau < 0.1 || fTau > 5.0) { - Warning("isTestAvailable", "On/Off model parameter 'tau' must be in the range [0.1,5.0]. Skipping test..."); - return false; - } - return true; - } - - bool testCode() override - { - - // names of tested variables must be the same in write / comparison modes - TString significanceString = TString::Format("thtc1_significance_%d_%d_%lf", fObsValueOn, fObsValueOff, fTau); - - // build workspace and model - auto w = std::make_unique("w"); - buildOnOffModel(*w); - ModelConfig *sbModel = static_cast(w->obj("S+B")); - ModelConfig *bModel = static_cast(w->obj("B")); - - // add observable values to data set and fix other parameters - w->var("n_on")->setVal(fObsValueOn); - w->var("n_off")->setVal(fObsValueOff); - w->var("tau")->setVal(fTau); - w->var("tau")->setConstant(); - w->data("data")->add(*sbModel->GetObservables()); - w->var("bkg")->setVal(fObsValueOff / fTau); - - // Make snapshots - w->var("sig")->setVal(fObsValueOn - fObsValueOff / fTau); - sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); - w->var("sig")->setVal(0.0); - bModel->SetSnapshot(*bModel->GetParametersOfInterest()); - - // Do hypothesis test with ProfileLikelihoodCalculator - if (_write == true) { - - auto plc = std::make_unique(*w->data("data"), *sbModel); - plc->SetNullParameters(*bModel->GetSnapshot()); - plc->SetAlternateParameters(*sbModel->GetSnapshot()); - regValue(plc->GetHypoTest()->Significance(), significanceString); - - } else { // Do hypothesis test with AsymptoticCalculator - - AsymptoticCalculator::SetPrintLevel(_verb); // disable superfluous messaging - auto atc = std::make_unique(*w->data("data"), *sbModel, *bModel); - atc->SetOneSidedDiscovery(true); - regValue(atc->GetHypoTest()->Significance(), significanceString); - } - - return true; - } -}; - -/////////////////////////////////////////////////////////////////////////////// -// -// HYPOTHESIS TEST CALCULATOR TEST - SIMULTANEOUS PDF MODEL -// -// This test evaluates the functionality of the HypoTestCalculator by -// calculating the significance of the signal on a simple Simultaneous Pdf -// model with two channels. Reference values and test values are both computed -// with the HypoTestCalculator. As such, this test can only confirm if the -// HypoTestCalculator has the same behaviour across different computing -// platforms or RooStats revisions. -// -// ModelConfig (explicit) : Simultaneous Model -// built in stressRooStats_models.cxx -// -// Input Parameters: -// calculatorType -> Frequentist, Hybrid or Asymptotic -// testStatType -> Profile Likelihood Ratio, Simple Likelihood Ratio, etc... -// -// 06/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestHypoTestCalculator2 : public RooUnitTest { -private: - ECalculatorType fCalculatorType; - ETestStatType fTestStatType; - -public: - TestHypoTestCalculator2(TFile *refFile, bool writeRef, Int_t verbose, ECalculatorType calculatorType = kAsymptotic, - ETestStatType testStatType = kProfileLROneSidedDiscovery) - : RooUnitTest(TString::Format("HypoTestCalculator Significance - Simultaneous Pdf - %s - %s", - kECalculatorTypeString[calculatorType], kETestStatTypeString[testStatType]), - refFile, writeRef, verbose), - fCalculatorType(calculatorType), - fTestStatType(testStatType){}; - - bool testCode() override - { - - // Build workspace and models - auto w = std::make_unique("w"); - buildSimultaneousModel(w.get()); - ModelConfig *sbModel = (ModelConfig *)w->obj("S+B"); - ModelConfig *bModel = (ModelConfig *)w->obj("B"); - - // set snapshots - sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); // value set in model - w->var("sig")->setVal(0); - bModel->SetSnapshot(*bModel->GetParametersOfInterest()); - - AsymptoticCalculator::SetPrintLevel(_verb); // is static (don;t care if we don't use it) - - std::unique_ptr calc{buildHypoTestCalculator(fCalculatorType, *w->data("data"), *bModel, *sbModel, 500, 50)}; - if (fCalculatorType == kAsymptotic) { - ((AsymptoticCalculator *)calc.get())->SetOneSidedDiscovery(true); - } - - // ToyMCSampler configuration - ToyMCSampler *tmcs = (ToyMCSampler *)calc->GetTestStatSampler(); - tmcs->SetTestStatistic(buildTestStatistic(fTestStatType, *bModel, *sbModel)); - tmcs->SetUseMultiGen(true); // speedup - - // Register result (test significance) - std::unique_ptr htr{calc->GetHypoTest()}; - regValue(htr->Significance(), TString::Format("thtc2_significance_%s_%s", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType])); - - // corresponding visual plots (in verbose mode) - from tutorials/roofit/roostats/StandardHypoTestDemo.C - if (_verb >= 1) { - if (fCalculatorType != kAsymptotic) { - TCanvas *c = new TCanvas("thtc2_canvas", "THTC2 Canvas"); - - c->cd(1); - HypoTestPlot *plot = new HypoTestPlot(*htr, 100); - plot->SetLogYaxis(true); - plot->Draw(); - - SamplingDistribution *altDist = htr->GetAltDistribution(); - HypoTestResult htExp("Expected result"); - htExp.Append(htr.get()); - // find quantiles in alt (S+B) distribution - double p[5]; - double q[5]; - for (Int_t i = 0; i < 5; ++i) { - double sig = -2 + i; - p[i] = ROOT::Math::normal_cdf(sig, 1); - } - std::vector values = altDist->GetSamplingDistribution(); - TMath::Quantiles(values.size(), 5, &values[0], q, p, false); - - for (Int_t i = 0; i < 5; ++i) { - htExp.SetTestStatisticData(q[i]); - double sig = -2 + i; - std::cout << "Expected p-value and significance at " << sig << " sigma = " << htExp.NullPValue() - << " significance " << htExp.Significance() << " sigma " << std::endl; - } - c->SaveAs(TString::Format("thtc2_scan_%s_%s.pdf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType])); - } else { - for (Int_t i = 0; i < 5; ++i) { - double sig = -2 + i; - double pval = - AsymptoticCalculator::GetExpectedPValues(htr->NullPValue(), htr->AlternatePValue(), -sig, false); - std::cout << "Expected p-value and significance at " << sig << " sigma = " << pval << " significance " - << ROOT::Math::normal_quantile_c(pval, 1) << " sigma " << std::endl; - } - } - } - - return true; - } -}; - -// -// END OF PART FOUR -// -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -//_____________________________________________________________________________ - -//_____________________________________________________________________________ -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -// -// PART FIVE: -// HYPOTHESIS TEST INVERTER UNIT TESTS -// - -/////////////////////////////////////////////////////////////////////////////// -// -// HYPOTESTINVERTER INTERVAL - POISSON PRODUCT MODEL -// -// Test the validity of the confidence interval computed by the HypoTestInverter -// on a complex Poisson model distribution. Reference values and test values -// are both computed with the HypoTestInverter. As such, this test can only -// confirm if the HypoTestInverter has the same behaviour across different -// computing platforms or RooStats revisions. -// -// ModelConfig (explicit) : Poisson Product Model -// built in stressRooStats_models.cxx -// -// Input Parameters: -// calculatorType -> Frequentist, Hybrid or Asymptotic -// testStatType -> Profile Likelihood Ratio, Simple Likelihood Ratio, etc... -// obsValueX -> observed value "x" when measuring sig + bkg1 -// obsValueY -> observed value "y" when measuring 2*sig*1.2^beta + bkg2 -// confidenceLevel -> Confidence Level of the interval we are calculating -// -// 04/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestHypoTestInverter1 : public RooUnitTest { -private: - ECalculatorType fCalculatorType; - ETestStatType fTestStatType; - Int_t fObsValueX; - Int_t fObsValueY; - double fConfidenceLevel; - -public: - TestHypoTestInverter1(TFile *refFile, bool writeRef, Int_t verbose, ECalculatorType calculatorType = kAsymptotic, - ETestStatType testStatType = kProfileLR, Int_t obsValueX = 15, Int_t obsValueY = 30, - double confidenceLevel = 2 * normal_cdf(1) - 1) - : RooUnitTest(TString::Format("HypoTestInverter Interval - Poisson Product Model - %s - %s", - kECalculatorTypeString[calculatorType], kETestStatTypeString[testStatType]), - refFile, writeRef, verbose), - fCalculatorType(calculatorType), - fTestStatType(testStatType), - fObsValueX(obsValueX), - fObsValueY(obsValueY), - fConfidenceLevel(confidenceLevel){}; - - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fObsValueX < 0 || fObsValueX > 30) { - Warning("isTestAvailable", "Observed value X=s+b must be in the range [0,30]. Skipping test..."); - return false; - } - if (fObsValueY < 0 || fObsValueY > 80) { - Warning("isTestAvailable", "Observed value Y=2*s*1.2^beta+b must be in the range [0,80]. Skipping test..."); - return false; - } - if (fConfidenceLevel <= 0.0 || fConfidenceLevel >= 1.0) { - Warning("isTestAvailable", "Confidence level must be in the range (0,1). Skipping test..."); - return false; - } - return true; - } - - // larger value test tolerance especially when using toys (difference of <~ 0.1 observed between using Minuit or - // Minuit2) - // (inherited default value is 1e-3) - double vtol() override { return (fCalculatorType == kAsymptotic) ? 0.01 : 0.1; } - - bool testCode() override - { - - // Create workspace and model - auto w = std::make_unique("w"); - buildPoissonProductModel(w.get()); - ModelConfig *sbModel = (ModelConfig *)w->obj("S+B"); - ModelConfig *bModel = (ModelConfig *)w->obj("B"); - - // add observed values to data set - w->var("x")->setVal(fObsValueX); - w->var("y")->setVal(fObsValueY); - w->data("data")->add(*sbModel->GetObservables()); - - std::unique_ptr initialVariables{sbModel->GetPdf()->getVariables()}; - w->saveSnapshot("initialVariables", *initialVariables); - - // set snapshots - w->var("sig")->setVal(fObsValueX - w->var("bkg1")->getValV()); - sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); - w->var("sig")->setVal(0); - bModel->SetSnapshot(*bModel->GetParametersOfInterest()); - - // build and configure HypoTestInverter - AsymptoticCalculator::SetPrintLevel(_verb); - HypoTestCalculatorGeneric *calc = - buildHypoTestCalculator(fCalculatorType, *w->data("data"), *sbModel, *bModel, 100, 1); - auto hti = std::make_unique(*calc, nullptr, 1.0 - fConfidenceLevel); - hti->SetTestStatistic(*buildTestStatistic(fTestStatType, *sbModel, *bModel)); - hti->SetVerbose(_verb); - - int nscanPoints = 10; - if (fCalculatorType == kAsymptotic) { - ((AsymptoticCalculator *)calc)->SetTwoSided(); - ((AsymptoticCalculator *)calc)->SetPrintLevel(_verb); - nscanPoints = 40; - } - - hti->SetFixedScan(nscanPoints, w->var("sig")->getMin(), w->var("sig")->getMax()); // significant speedup - - // ToyMCSampler configuration - ToyMCSampler *tmcs = (ToyMCSampler *)hti->GetHypoTestCalculator()->GetTestStatSampler(); - tmcs->SetNEventsPerToy(1); // needed because we don't have an extended pdf - tmcs->SetUseMultiGen(true); // speedup - - std::unique_ptr interval{hti->GetInterval()}; - regValue(interval->LowerLimit(), - TString::Format("thti1_lower_limit_sig_%s_%s_%d_%d_%lf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fObsValueY, fConfidenceLevel)); - regValue(interval->UpperLimit(), - TString::Format("thti1_upper_limit_sig_%s_%s_%d_%d_%lf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fObsValueY, fConfidenceLevel)); - - if (_verb >= 1) { - HypoTestInverterPlot *plot = new HypoTestInverterPlot("thti1_scan", "Two-Sided Scan", interval.get()); - TCanvas *c1 = new TCanvas("thti1_canvas", "THTI Canvas"); - c1->SetLogy(false); - plot->Draw("2CL CLB"); - c1->SaveAs(TString::Format("thti1_scan_%s_%s_%d_%d_%lf.pdf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fObsValueY, fConfidenceLevel)); - - if (_verb == 2) { - const int n = interval->ArraySize(); - if (n > 0 && interval->GetResult(0)->GetNullDistribution()) { - TCanvas *c2 = new TCanvas("thti1_teststat_dist", "HTI Test Statistic Distributions", 2); - if (n > 1) { - int ny = TMath::CeilNint(sqrt((double)n)); - int nx = TMath::CeilNint(double(n) / ny); - c2->Divide(nx, ny); - } - for (int i = 0; i < n; ++i) { - if (n > 1) - c2->cd(i + 1); - SamplingDistPlot *pl = plot->MakeTestStatPlot(i); - if (pl == nullptr) - return true; - pl->SetLogYaxis(true); - pl->Draw(); - } - c2->SaveAs(TString::Format("thti1_teststat_distrib_%s_%s_%d_%d_%lf.pdf", - kECalculatorTypeString[fCalculatorType], kETestStatTypeString[fTestStatType], - fObsValueX, fObsValueY, fConfidenceLevel)); - } - } - } - - // in case of debug write the workspace in a file - if (_verb > 1) { - w->loadSnapshot("initialVariables"); - w->writeToFile(TString::Format("stressRooStats_PoissonProductModel_%d_%d.root", fObsValueX, fObsValueY)); - } - - return true; - } -}; - -/////////////////////////////////////////////////////////////////////////////// -// -// HYPOTESTINVERTER UPPER LIMIT - SIGNAL + BACKGROUND + EFFICIENCY MODEL -// -// Test the validity of the upper limit computed by the HypoTestInverter -// on a complex model distribution with signal, background and efficiency. -// Reference values and test values are both computed with the HypoTestInverter. -// As such, this test can only confirm if the HypoTestInverter has the same -// behaviour across different computing platforms or RooStats revisions. -// -// ModelConfig (explicit) : Poisson Signal + Background + Efficiency -// built in stressRooStats_models.cxx -// -/// Input Parameters: -// calculatorType -> Frequentist, Hybrid or Asymptotic -// testStatType -> Profile Likelihood Ratio, Simple Likelihood Ratio, etc... -// obsValueX -> observed value "x" when measuring sig * eff + bkg -// confidenceLevel -> Confidence Level of the upper limit we are calculating -// -// 04/2012 - Ioan Gabriel Bucur -// -/////////////////////////////////////////////////////////////////////////////// - -class TestHypoTestInverter2 : public RooUnitTest { -private: - ECalculatorType fCalculatorType; - ETestStatType fTestStatType; - Int_t fObsValueX; - double fConfidenceLevel; - -public: - TestHypoTestInverter2(TFile *refFile, bool writeRef, Int_t verbose, ECalculatorType calculatorType = kAsymptotic, - ETestStatType testStatType = kProfileLROneSided, Int_t obsValueX = 10, - double confidenceLevel = 2 * normal_cdf(1) - 1) - : RooUnitTest(TString::Format("HypoTestInverter Upper Limit - Poisson Efficiency Model - %s - %s", - kECalculatorTypeString[calculatorType], kETestStatTypeString[testStatType]), - refFile, writeRef, verbose), - fCalculatorType(calculatorType), - fTestStatType(testStatType), - fObsValueX(obsValueX), - fConfidenceLevel(confidenceLevel){}; - - // larger value test tolerance especially when using toys (difference of <~ 0.1 observed between using Minuit or - // Minuit2) - // (inherited default value is 1e-3) - double vtol() override { return (fCalculatorType == kAsymptotic) ? 0.02 : 0.1; } - // Basic checks for the parameters passed to the test - // In case of invalid parameters, a warning is printed and the test is skipped - bool isTestAvailable() override - { - if (fObsValueX < 0 || fObsValueX > 50) { - Warning("isTestAvailable", "Observed value X=s*e+b must be in the range [0,70]. Skipping test..."); - return false; - } - if (fConfidenceLevel <= 0.0 || fConfidenceLevel >= 1.0) { - Warning("isTestAvailable", "Confidence level must be in the range (0,1). Skipping test..."); - return false; - } - return true; - } - - bool testCode() override - { - - // Create workspace and model - auto w = std::make_unique("w"); - buildPoissonEfficiencyModel(*w); - ModelConfig *sbModel = (ModelConfig *)w->obj("S+B"); - ModelConfig *bModel = (ModelConfig *)w->obj("B"); - - // add observed values to data set - w->var("x")->setVal(fObsValueX); - w->data("data")->add(*sbModel->GetObservables()); - - std::unique_ptr initialVariables{sbModel->GetPdf()->getVariables()}; - w->saveSnapshot("initialVariables", *initialVariables); - - // set snapshots - sbModel->SetSnapshot(*sbModel->GetParametersOfInterest()); - w->var("sig")->setVal(0); - bModel->SetSnapshot(*bModel->GetParametersOfInterest()); - - // calculate upper limit with HypoTestInverter - AsymptoticCalculator::SetPrintLevel(_verb); - HypoTestCalculatorGeneric *calc = - buildHypoTestCalculator(fCalculatorType, *w->data("data"), *sbModel, *bModel, 100, 100); - auto hti = std::make_unique(*calc, nullptr, 1.0 - fConfidenceLevel); - hti->SetTestStatistic(*buildTestStatistic(fTestStatType, *sbModel, *bModel)); - hti->SetVerbose(_verb); - - int npoints = 10; - if (fCalculatorType == kAsymptotic) { - ((AsymptoticCalculator *)calc)->SetOneSided(true); - ((AsymptoticCalculator *)calc)->SetPrintLevel(_verb); - npoints = 40; - } - - hti->SetFixedScan(npoints, w->var("sig")->getMin(), w->var("sig")->getMax()); // significant speedup - - // needed because we have no extended pdf and the ToyMC Sampler evaluation returns an error - ToyMCSampler *tmcs = (ToyMCSampler *)hti->GetHypoTestCalculator()->GetTestStatSampler(); - tmcs->SetNEventsPerToy(1); - tmcs->SetUseMultiGen(true); // make ToyMCSampler faster - - // calculate interval and extract observed upper limit and expected upper limit (+- sigma) - std::unique_ptr interval{hti->GetInterval()}; - regValue(interval->UpperLimit(), - TString::Format("thti2_upper_limit_sig_%s_%s_%d_%lf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fConfidenceLevel)); - regValue(interval->GetExpectedUpperLimit(0), - TString::Format("thti2_exp_upper_limit_sig_%s_%s_%d_%lf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fConfidenceLevel)); - regValue(interval->GetExpectedUpperLimit(-2), - TString::Format("thti2_exp_upper_limit_-2_sig_%s_%s_%d_%lf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fConfidenceLevel)); - regValue(interval->GetExpectedUpperLimit(-1), - TString::Format("thti2_exp_upper_limit_-1_sig_%s_%s_%d_%lf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fConfidenceLevel)); - regValue(interval->GetExpectedUpperLimit(1), - TString::Format("thti2_exp_upper_limit_+1_sig_%s_%s_%d_%lf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fConfidenceLevel)); - regValue(interval->GetExpectedUpperLimit(2), - TString::Format("thti2_exp_upper_limit_+2_sig_%s_%s_%d_%lf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fConfidenceLevel)); - - if (_verb >= 1) { - HypoTestInverterPlot *plot = new HypoTestInverterPlot("thti2_scan", "HTI Upper Limit Scan", interval.get()); - TCanvas *c1 = new TCanvas("HypoTestInverter Scan"); - c1->SetLogy(false); - plot->Draw("2CL CLB"); - c1->SaveAs(TString::Format("thti2_scan_%s_%s_%d_%lf.pdf", kECalculatorTypeString[fCalculatorType], - kETestStatTypeString[fTestStatType], fObsValueX, fConfidenceLevel)); - - if (_verb == 2) { - const int n = interval->ArraySize(); - if (n > 0 && interval->GetResult(0)->GetNullDistribution()) { - TCanvas *c2 = new TCanvas("thti2_teststat_dist", "HTI Test Statistic Distributions", 2); - if (n > 1) { - int ny = TMath::CeilNint(sqrt((double)n)); - int nx = TMath::CeilNint(double(n) / ny); - c2->Divide(nx, ny); - } - for (int i = 0; i < n; ++i) { - if (n > 1) - c2->cd(i + 1); - SamplingDistPlot *pl = plot->MakeTestStatPlot(i); - if (pl == nullptr) - return true; - pl->SetLogYaxis(true); - pl->Draw(); - } - c2->SaveAs(TString::Format("thti2_teststat_distrib_%s_%s_%d_%lf.pdf", - kECalculatorTypeString[fCalculatorType], kETestStatTypeString[fTestStatType], - fObsValueX, fConfidenceLevel)); - } - } - } - - if (_verb > 1) { - w->loadSnapshot("initialVariables"); - w->writeToFile("stressRooStats_PoissonEfficiencyModel.root"); - } - - return true; - } -}; - -// -// END OF PART FIVE -// -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// -//_____________________________________________________________________________ - -// Other tests currently not included in any suite - -class TestHypoTestCalculator : public RooUnitTest { -public: - TestHypoTestCalculator(TFile *refFile, bool writeRef, Int_t verbose) - : RooUnitTest("HypoTestCalculator - On / Off Problem", refFile, writeRef, verbose){}; - - bool testCode() override - { - - const Int_t xValue = 150; - const Int_t yValue = 100; - const double tauValue = 1.0; - - if (_write == true) { - - // register analytical Z_Bi value - double Z_Bi = NumberCountingUtils::BinomialWithTauObsZ(xValue, yValue, tauValue); - regValue(Z_Bi, "thtc_significance_hybrid"); - - } else { - - // Make model for prototype on/off problem - // Pois(x | s+b) * Pois(y | tau b ) - auto w = std::make_unique("w"); - w->factory( - TString::Format("Poisson::on_pdf(x[%d,0,500],sum::splusb(sig[0,0,100],bkg[100,0,300]))", xValue).Data()); - w->factory(TString::Format("Poisson::off_pdf(y[%d,0,500],prod::taub(tau[%lf],bkg))", yValue, tauValue).Data()); - w->factory("PROD::prod_pdf(on_pdf, off_pdf)"); - - w->var("x")->setVal(xValue); - w->var("y")->setVal(yValue); - w->var("y")->setConstant(); - w->var("tau")->setVal(tauValue); - - // construct the Bayesian-averaged model (eg. a projection pdf) - // p'(x|s) = \int db p(x|s+b) * [ p(y|b) * prior(b) ] - w->factory("Uniform::prior(bkg)"); - w->factory("PROJ::averagedModel(PROD::foo(on_pdf|bkg,off_pdf,prior),bkg)"); - - // define sets of variables obs={x} and poi={sig} - // x is the only observable in the main measurement and y is treated as a separate measurement, - // which is used to produce the prior that will be used in the calculation to randomize the nuisance parameters - w->defineSet("obs", "x"); - w->defineSet("poi", "sig"); - - // Add observable value to a data set - RooDataSet data{"data", "data", *w->set("obs")}; - data.add(*w->set("obs")); - - // Build S+B and B models - ModelConfig *sbModel = new ModelConfig("SB_ModelConfig", w.get()); - sbModel->SetPdf(*w->pdf("prod_pdf")); - sbModel->SetObservables(*w->set("obs")); - sbModel->SetParametersOfInterest(*w->set("poi")); - w->var("sig")->setVal(xValue - yValue / tauValue); // important ! - sbModel->SetSnapshot(*w->set("poi")); - - ModelConfig *bModel = new ModelConfig("B_ModelConfig", w.get()); - bModel->SetPdf(*w->pdf("prod_pdf")); - bModel->SetObservables(*w->set("obs")); - bModel->SetParametersOfInterest(*w->set("poi")); - w->var("sig")->setVal(0.0); // important ! - bModel->SetSnapshot(*w->set("poi")); - - // alternate priors - w->factory("Gaussian::gauss_prior(bkg, y, expr::sqrty('sqrt(y)', y))"); - w->factory("Lognormal::lognorm_prior(bkg, y, expr::kappa('1+1./sqrt(y)',y))"); - - // build test statistic - SimpleLikelihoodRatioTestStat *slrts = - new SimpleLikelihoodRatioTestStat(*bModel->GetPdf(), *sbModel->GetPdf()); - slrts->SetNullParameters(*bModel->GetSnapshot()); - slrts->SetAltParameters(*sbModel->GetSnapshot()); - slrts->SetAlwaysReuseNLL(true); - - RatioOfProfiledLikelihoodsTestStat *roplts = - new RatioOfProfiledLikelihoodsTestStat(*bModel->GetPdf(), *sbModel->GetPdf()); - roplts->SetAlwaysReuseNLL(true); - - ProfileLikelihoodTestStat *pllts = new ProfileLikelihoodTestStat(*bModel->GetPdf()); - pllts->SetAlwaysReuseNLL(true); - - MaxLikelihoodEstimateTestStat *mlets = new MaxLikelihoodEstimateTestStat( - *sbModel->GetPdf(), *((RooRealVar *)sbModel->GetParametersOfInterest()->first())); - - NumEventsTestStat *nevts = new NumEventsTestStat(*sbModel->GetPdf()); - - auto htc = std::make_unique(data, *sbModel, *bModel); - ToyMCSampler *tmcs = (ToyMCSampler *)htc->GetTestStatSampler(); - tmcs->SetNEventsPerToy(1); - htc->SetToys(5000, 1000); - htc->ForcePriorNuisanceAlt(*w->pdf("off_pdf")); - htc->ForcePriorNuisanceNull(*w->pdf("off_pdf")); - - tmcs->SetTestStatistic(pllts); - std::unique_ptr htr{htc->GetHypoTest()}; - htr->Print(); - std::cout << "PLLTS " << htr->Significance() << std::endl; - tmcs->SetTestStatistic(mlets); - htr = std::unique_ptr{htc->GetHypoTest()}; - htr->Print(); - std::cout << "MLETS " << htr->Significance() << std::endl; - tmcs->SetTestStatistic(nevts); - htr = std::unique_ptr{htc->GetHypoTest()}; - htr->Print(); - std::cout << "NEVTS " << htr->Significance() << std::endl; - tmcs->SetTestStatistic(slrts); - htr = std::unique_ptr{htc->GetHypoTest()}; - htr->Print(); - std::cout << "SLRTS " << htr->Significance() << std::endl; - tmcs->SetTestStatistic(roplts); - htr = std::unique_ptr{htc->GetHypoTest()}; - htr->Print(); - std::cout << "ROPLTS " << htr->Significance() << std::endl; - - regValue(htr->Significance(), "thtc_significance_hybrid"); - - if (_verb > 1) - w->writeToFile("stressRooStats_OnOffModel.root"); - } - - return true; - } -}; - -static HypoTestCalculatorGeneric *buildHypoTestCalculator(const ECalculatorType calculatorType, RooAbsData &data, - const ModelConfig &nullModel, const ModelConfig &altModel, - const UInt_t toysNull, const UInt_t toysAlt) -{ - HypoTestCalculatorGeneric *calc = nullptr; - - if (calculatorType == kAsymptotic) { - AsymptoticCalculator *ac = new AsymptoticCalculator(data, altModel, nullModel); - calc = ac; - } else if (calculatorType == kFrequentist) { - FrequentistCalculator *fc = new FrequentistCalculator(data, altModel, nullModel); - // set toys for speedup - fc->SetToys(toysNull, toysAlt); - calc = fc; - } else { // kHybrid - HybridCalculator *hc = new HybridCalculator(data, altModel, nullModel); - // set toys for speedup - hc->SetToys(toysNull, toysAlt); - calc = hc; - } - - assert(calc != NULL); // sanity check - should never happen - - return calc; -} - -static TestStatistic * -buildTestStatistic(const ETestStatType testStatType, const ModelConfig &nullModel, const ModelConfig &altModel) -{ - - TestStatistic *testStat = nullptr; - - if (testStatType == kSimpleLR) { - auto *slrts = new SimpleLikelihoodRatioTestStat(*nullModel.GetPdf(), *altModel.GetPdf()); - // TODO - different for HypoTestInverter and HypoTestCalculator - RooArgSet nullParams(*nullModel.GetSnapshot()); - if (nullModel.GetNuisanceParameters()) - nullParams.add(*nullModel.GetNuisanceParameters()); - if (nullModel.GetSnapshot()) - slrts->SetNullParameters(nullParams); - RooArgSet altParams(*altModel.GetSnapshot()); - if (altModel.GetNuisanceParameters()) - altParams.add(*altModel.GetNuisanceParameters()); - if (altModel.GetSnapshot()) - slrts->SetAltParameters(altParams); - slrts->SetAlwaysReuseNLL(true); - testStat = slrts; - } else if (testStatType == kRatioLR) { - auto *roplts = - new RatioOfProfiledLikelihoodsTestStat(*nullModel.GetPdf(), *altModel.GetPdf(), altModel.GetSnapshot()); - roplts->SetSubtractMLE(false); - roplts->SetAlwaysReuseNLL(true); - testStat = roplts; - } else if (testStatType == kMLE) { - auto *mlets = new MaxLikelihoodEstimateTestStat(*nullModel.GetPdf(), - *((RooRealVar *)nullModel.GetParametersOfInterest()->first())); - testStat = mlets; - } else if (testStatType == kNObs) { - NumEventsTestStat *nevtts = new NumEventsTestStat(*nullModel.GetPdf()); - testStat = nevtts; - } else { // kProfileLR, kProfileLROneSided and kProfileLRSigned - auto *plts = new ProfileLikelihoodTestStat(*nullModel.GetPdf()); - if (testStatType == kProfileLROneSided) { - plts->SetOneSided(true); - } else if (testStatType == kProfileLROneSidedDiscovery) { - plts->SetOneSidedDiscovery(true); - } else if (testStatType == kProfileLRSigned) { - plts->SetSigned(true); - } - plts->SetAlwaysReuseNLL(true); - testStat = plts; - } - - assert(testStat != nullptr); // sanity check - should never happen - - return testStat; -}