Skip to content
Snippets Groups Projects
Commit 8d8688ef authored by Lorenzo Moneta's avatar Lorenzo Moneta
Browse files

from Gabriel: new version with tests on HypoTestInverter 

git-svn-id: http://root.cern.ch/svn/root/trunk@44165 27541ba8-7e3a-0410-8455-c3a389f83636
parent 36ae1b1e
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
test/stressRooStats.cxx
+ 45
− 22
View file @ 8d8688ef
// @(#)root/roofitcore:$name: $:$id$
// Authors: Wouter Verkerke November 2007
// C/C++ headers
#include <string>
#include <list>
#include <iostream>
#include <cmath>
// ROOT headers
#include "TWebFile.h"
#include "TSystem.h"
#include "TString.h"
......@@ -13,6 +20,8 @@
#include "TH1.h"
#include "TF1.h"
#include "TBenchmark.h"
// RooFit headers
#include "RooGlobalFunc.h"
#include "RooNumIntConfig.h"
#include "RooMsgService.h"
......@@ -27,12 +36,8 @@
#include "RooRandom.h"
#include "RooTrace.h"
#include "RooMath.h"
#include <string>
#include <list>
#include <iostream>
#include <math.h>
// Tests file
#include "stressRooStats_tests.cxx"
using namespace std ;
......@@ -55,7 +60,6 @@ void StatusPrint(const Int_t id, const TString &title, const Int_t status, const
cout << left << setw(lineWidth) << setfill('.') << header << " " << (status > 0 ? "OK" : (status < 0 ? "SKIPPED" : "FAILED")) << endl;
}
//______________________________________________________________________________
Int_t stressRooStats(const char* refFile, Bool_t writeRef, Int_t verbose, Int_t oneTest, Bool_t dryRun, Bool_t doDump, Bool_t doTreeStore)
{
......@@ -125,8 +129,8 @@ Int_t stressRooStats(const char* refFile, Bool_t writeRef, Int_t verbose, Int_t
// TEST PLC CONFINT PRODUCT POISSON : Observed value range is [0,40] for x=s+b and [0,120] 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, 25, 40));
testList.push_back(new TestProfileLikelihoodCalculator3(fref, writeRef, verbose, 20, 25, 2 * ROOT::Math::normal_cdf(2) - 1));
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));
// TEST PLC HYPOTEST ON/OFF MODEL
testList.push_back(new TestProfileLikelihoodCalculator4(fref, writeRef, verbose));
......@@ -142,28 +146,47 @@ Int_t stressRooStats(const char* refFile, Bool_t writeRef, Int_t verbose, Int_t
// TEST BC CONFINT CENTRAL PRODUCT POISSON : Observed value range is [0,40] for x=s+b and [0,120] 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, 25, 40));
testList.push_back(new TestBayesianCalculator3(fref, writeRef, verbose, 20, 25, 2 * ROOT::Math::normal_cdf(2) - 1));
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));
// TEST MCMCC CONFINT PRODUCT POISSON : Observed value range is [0,40] for x=s+b and [0,120] 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, 25, 40));
testList.push_back(new TestMCMCCalculator(fref, writeRef, verbose, 20, 25, 2 * ROOT::Math::normal_cdf(2) - 1));
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));
// TEST HYPOTEST INVERTER PRODUCT POISSON : Observed value range is [0,40] for x=s+b and [0,120] for y=2*s*1.2^beta
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kProfileLR, 10, 30));
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kProfileLR, 20, 25));
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kProfileLR, 15, 20, 2 * normal_cdf(2) - 1));
// TEST PLC CONFINT PRODUCT POISSON : Observed value range is [0,40] for x=s+b and [0,120] 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));
// testList.push_back(new TestBasic101(fref, writeRef, verbose));
// testList.push_back(new TestHypoTestCalculator(fref, writeRef, verbose));
// testList.push_back(new TestHypoTestCalculator2(fref, writeRef, verbose));
// testList.push_back(new TestHypoTestCalculator3(fref, writeRef, verbose));
// testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kAsymptotic, kSimpleLR));
// testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kHybrid, kSimpleLR));
// testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kSimpleLR));
// testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kRatioLR));
// testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kProfileLR));
// testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kProfileLROneSided));
// testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kProfileLRSigned));
// testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kMLE));
// testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kNObs));
/* testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kProfileLR, 10, 30));
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kProfileLR, 20, 25));
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kProfileLR, 15, 20, 2 * normal_cdf(2) - 1));
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kSimpleLR, 10, 30));
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kSimpleLR, 20, 25));
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kSimpleLR, 15, 20, 2 * normal_cdf(2) - 1));
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kRatioLR, 10, 30));
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kRatioLR, 20, 25));
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kRatioLR, 15, 20, 2 * normal_cdf(2) - 1));
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kMLE, 10, 30));
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kMLE, 20, 25));
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kMLE, 15, 20, 2 * normal_cdf(2) - 1));
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kNEvts, 10, 30));
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kRatioLR, 20, 25));
testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kFrequentist, kRatioLR, 15, 20, 2 * normal_cdf(2) - 1));
*/
// testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kHybrid, kSimpleLR));
// testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kHybrid, kRatioLR));
// testList.push_back(new TestHypoTestInverter1(fref, writeRef, verbose, HypoTestInverter::kHybrid, kProfileLR));
......
test/stressRooStats_models.cxx
+ 48
− 2
View file @ 8d8688ef
......@@ -2,6 +2,7 @@
#include "RooWorkspace.h"
#include "RooRealVar.h"
#include "RooDataSet.h"
// RooStats headers
#include "RooStats/ModelConfig.h"
......@@ -49,8 +50,6 @@ void buildPoissonProductModel(RooWorkspace *w)
// create background model configuration
ModelConfig *bModel = new ModelConfig(*sbModel);
bModel->SetName("B");
w->var("sig")->setVal(0);
bModel->SetSnapshot(*w->set("poi"));
w->import(*sbModel);
w->import(*bModel);
......@@ -100,3 +99,50 @@ void buildOnOffModel(RooWorkspace *w)
w->import(*bModel);
}
void createPoissonEfficiencyModel(RooWorkspace *w) {
// build models
w->factory("Gaussian::constrb(b0[-5,5], b1[-5,5], 1)");
w->factory("Gaussian::constre(e0[-5,5], e1[-5,5], 1)");
w->factory("expr::bkg('5 * pow(1.3, b1)', b1)"); // background model
w->factory("expr::eff('0.5 * pow(1.2, e1)', e1)"); // efficiency model
w->factory("expr::splusb('eff * sig + bkg', eff, bkg, sig[0,20])");
w->factory("Poisson::sb_poiss(x[0,40], splusb)");
w->factory("Poisson::b_poiss(x, bkg)");
w->factory("PROD::sb_pdf(sb_poiss, constrb, constre)");
w->factory("PROD::b_pdf(b_poiss, constrb)");
w->factory("PROD::priorbkg(constr1, constr2)");
w->var("b0")->setConstant(kTRUE);
w->var("e0")->setConstant(kTRUE);
// build argument sets
w->defineSet("obs", "x");
w->defineSet("poi", "sig");
w->defineSet("nuis", "b1,e1");
w->defineSet("globObs", "b0,e0");
// define data set and import it into workspace
RooDataSet *data = new RooDataSet("data", "data", *w->set("obs"));
w->import(*data);
// create model configuration
ModelConfig *sbModel = new ModelConfig("S+B", w);
sbModel->SetObservables(*w->set("obs"));
sbModel->SetParametersOfInterest(*w->set("poi"));
sbModel->SetNuisanceParameters(*w->set("nuis"));
sbModel->SetPdf("sb_pdf");
//sbModel->SetPriorPdf("prior");
sbModel->SetSnapshot(*w->set("poi"));
sbModel->SetGlobalObservables(*w->set("globObs"));
ModelConfig *bModel = new ModelConfig(*sbModel);
bModel->SetName("B");
bModel->SetPdf("b_pdf");
w->import(*sbModel);
w->import(*bModel);
}
test/stressRooStats_tests.cxx
+ 228
− 216
View file @ 8d8688ef
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment