// Example showing confidence intervals with four techniques. /* IntervalExamples Author Kyle Cranmer date Sep. 2010 An example that shows confidence intervals with four techniques. The model is a Normal Gaussian G(x|mu,sigma) with 100 samples of x. The answer is known analytically, so this is a good example to validate the RooStats tools. expected interval is [-0.162917, 0.229075] plc interval is [-0.162917, 0.229075] fc interval is [-0.17 , 0.23] // stepsize is 0.01 bc interval is [-0.162918, 0.229076] mcmc interval is [-0.166999, 0.230224] */ #include "RooStats/ConfInterval.h" #include "RooStats/PointSetInterval.h" #include "RooStats/ConfidenceBelt.h" #include "RooStats/FeldmanCousins.h" #include "RooStats/ProfileLikelihoodCalculator.h" #include "RooStats/MCMCCalculator.h" #include "RooStats/BayesianCalculator.h" #include "RooStats/MCMCIntervalPlot.h" #include "RooStats/LikelihoodIntervalPlot.h" #include "RooStats/ProofConfig.h" #include "RooStats/ToyMCSampler.h" #include "RooRandom.h" #include "RooDataSet.h" #include "RooRealVar.h" #include "RooConstVar.h" #include "RooAddition.h" #include "RooDataHist.h" #include "RooPoisson.h" #include "RooPlot.h" #include "TCanvas.h" #include "TTree.h" #include "TStyle.h" #include "TMath.h" #include"Math/DistFunc.h" #include "TH1F.h" #include "TMarker.h" #include "TStopwatch.h" #include <iostream> // use this order for safety on library loading using namespace RooFit ; using namespace RooStats ; void IntervalExamples() { // Time this macro TStopwatch t; t.Start(); // set RooFit random seed for reproducible results RooRandom::randomGenerator()->SetSeed(3001); // make a simple model via the workspace factory RooWorkspace* wspace = new RooWorkspace(); wspace->factory("Gaussian::normal(x[-10,10],mu[-1,1],sigma[1])"); wspace->defineSet("poi","mu"); wspace->defineSet("obs","x"); // specify components of model for statistical tools ModelConfig* modelConfig = new ModelConfig("Example G(x|mu,1)"); modelConfig->SetWorkspace(*wspace); modelConfig->SetPdf( *wspace->pdf("normal") ); modelConfig->SetParametersOfInterest( *wspace->set("poi") ); modelConfig->SetObservables( *wspace->set("obs") ); // create a toy dataset RooDataSet* data = wspace->pdf("normal")->generate(*wspace->set("obs"),100); data->Print(); // for convenience later on RooRealVar* x = wspace->var("x"); RooRealVar* mu = wspace->var("mu"); // set confidence level double confidenceLevel = 0.95; // example use profile likelihood calculator ProfileLikelihoodCalculator plc(*data, *modelConfig); plc.SetConfidenceLevel( confidenceLevel); LikelihoodInterval* plInt = plc.GetInterval(); // example use of Feldman-Cousins FeldmanCousins fc(*data, *modelConfig); fc.SetConfidenceLevel( confidenceLevel); fc.SetNBins(100); // number of points to test per parameter fc.UseAdaptiveSampling(true); // make it go faster // Here, we consider only ensembles with 100 events // The PDF could be extended and this could be removed fc.FluctuateNumDataEntries(false); // Proof // ProofConfig pc(*wspace, 4, "workers=4", kFALSE); // proof-lite //ProofConfig pc(w, 8, "localhost"); // proof cluster at "localhost" // ToyMCSampler* toymcsampler = (ToyMCSampler*) fc.GetTestStatSampler(); // toymcsampler->SetProofConfig(&pc); // enable proof PointSetInterval* interval = (PointSetInterval*) fc.GetInterval(); // example use of BayesianCalculator // now we also need to specify a prior in the ModelConfig wspace->factory("Uniform::prior(mu)"); modelConfig->SetPriorPdf(*wspace->pdf("prior")); // example usage of BayesianCalculator BayesianCalculator bc(*data, *modelConfig); bc.SetConfidenceLevel( confidenceLevel); SimpleInterval* bcInt = bc.GetInterval(); // example use of MCMCInterval MCMCCalculator mc(*data, *modelConfig); mc.SetConfidenceLevel( confidenceLevel); // special options mc.SetNumBins(200); // bins used internally for representing posterior mc.SetNumBurnInSteps(500); // first N steps to be ignored as burn-in mc.SetNumIters(100000); // how long to run chain mc.SetLeftSideTailFraction(0.5); // for central interval MCMCInterval* mcInt = mc.GetInterval(); // for this example we know the expected intervals double expectedLL = data->mean(*x) + ROOT::Math::normal_quantile( (1-confidenceLevel)/2,1) / sqrt(data->numEntries()); double expectedUL = data->mean(*x) + ROOT::Math::normal_quantile_c((1-confidenceLevel)/2,1) / sqrt(data->numEntries()) ; // Use the intervals std::cout << "expected interval is [" << expectedLL << ", " << expectedUL << "]" << endl; cout << "plc interval is [" << plInt->LowerLimit(*mu) << ", " << plInt->UpperLimit(*mu) << "]" << endl; std::cout << "fc interval is ["<< interval->LowerLimit(*mu) << " , " << interval->UpperLimit(*mu) << "]" << endl; cout << "bc interval is [" << bcInt->LowerLimit() << ", " << bcInt->UpperLimit() << "]" << endl; cout << "mc interval is [" << mcInt->LowerLimit(*mu) << ", " << mcInt->UpperLimit(*mu) << "]" << endl; mu->setVal(0); cout << "is mu=0 in the interval? " << plInt->IsInInterval(RooArgSet(*mu)) << endl; // make a reasonable style gStyle->SetCanvasColor(0); gStyle->SetCanvasBorderMode(0); gStyle->SetPadBorderMode(0); gStyle->SetPadColor(0); gStyle->SetCanvasColor(0); gStyle->SetTitleFillColor(0); gStyle->SetFillColor(0); gStyle->SetFrameFillColor(0); gStyle->SetStatColor(0); // some plots TCanvas* canvas = new TCanvas("canvas"); canvas->Divide(2,2); // plot the data canvas->cd(1); RooPlot* frame = x->frame(); data->plotOn(frame); data->statOn(frame); frame->Draw(); // plot the profile likeihood canvas->cd(2); LikelihoodIntervalPlot plot(plInt); plot.Draw(); // plot the MCMC interval canvas->cd(3); MCMCIntervalPlot* mcPlot = new MCMCIntervalPlot(*mcInt); mcPlot->SetLineColor(kGreen); mcPlot->SetLineWidth(2); mcPlot->Draw(); canvas->cd(4); RooPlot * bcPlot = bc.GetPosteriorPlot(); bcPlot->Draw(); canvas->Update(); t.Stop(); t.Print(); } IntervalExamples.C:1 IntervalExamples.C:2 IntervalExamples.C:3 IntervalExamples.C:4 IntervalExamples.C:5 IntervalExamples.C:6 IntervalExamples.C:7 IntervalExamples.C:8 IntervalExamples.C:9 IntervalExamples.C:10 IntervalExamples.C:11 IntervalExamples.C:12 IntervalExamples.C:13 IntervalExamples.C:14 IntervalExamples.C:15 IntervalExamples.C:16 IntervalExamples.C:17 IntervalExamples.C:18 IntervalExamples.C:19 IntervalExamples.C:20 IntervalExamples.C:21 IntervalExamples.C:22 IntervalExamples.C:23 IntervalExamples.C:24 IntervalExamples.C:25 IntervalExamples.C:26 IntervalExamples.C:27 IntervalExamples.C:28 IntervalExamples.C:29 IntervalExamples.C:30 IntervalExamples.C:31 IntervalExamples.C:32 IntervalExamples.C:33 IntervalExamples.C:34 IntervalExamples.C:35 IntervalExamples.C:36 IntervalExamples.C:37 IntervalExamples.C:38 IntervalExamples.C:39 IntervalExamples.C:40 IntervalExamples.C:41 IntervalExamples.C:42 IntervalExamples.C:43 IntervalExamples.C:44 IntervalExamples.C:45 IntervalExamples.C:46 IntervalExamples.C:47 IntervalExamples.C:48 IntervalExamples.C:49 IntervalExamples.C:50 IntervalExamples.C:51 IntervalExamples.C:52 IntervalExamples.C:53 IntervalExamples.C:54 IntervalExamples.C:55 IntervalExamples.C:56 IntervalExamples.C:57 IntervalExamples.C:58 IntervalExamples.C:59 IntervalExamples.C:60 IntervalExamples.C:61 IntervalExamples.C:62 IntervalExamples.C:63 IntervalExamples.C:64 IntervalExamples.C:65 IntervalExamples.C:66 IntervalExamples.C:67 IntervalExamples.C:68 IntervalExamples.C:69 IntervalExamples.C:70 IntervalExamples.C:71 IntervalExamples.C:72 IntervalExamples.C:73 IntervalExamples.C:74 IntervalExamples.C:75 IntervalExamples.C:76 IntervalExamples.C:77 IntervalExamples.C:78 IntervalExamples.C:79 IntervalExamples.C:80 IntervalExamples.C:81 IntervalExamples.C:82 IntervalExamples.C:83 IntervalExamples.C:84 IntervalExamples.C:85 IntervalExamples.C:86 IntervalExamples.C:87 IntervalExamples.C:88 IntervalExamples.C:89 IntervalExamples.C:90 IntervalExamples.C:91 IntervalExamples.C:92 IntervalExamples.C:93 IntervalExamples.C:94 IntervalExamples.C:95 IntervalExamples.C:96 IntervalExamples.C:97 IntervalExamples.C:98 IntervalExamples.C:99 IntervalExamples.C:100 IntervalExamples.C:101 IntervalExamples.C:102 IntervalExamples.C:103 IntervalExamples.C:104 IntervalExamples.C:105 IntervalExamples.C:106 IntervalExamples.C:107 IntervalExamples.C:108 IntervalExamples.C:109 IntervalExamples.C:110 IntervalExamples.C:111 IntervalExamples.C:112 IntervalExamples.C:113 IntervalExamples.C:114 IntervalExamples.C:115 IntervalExamples.C:116 IntervalExamples.C:117 IntervalExamples.C:118 IntervalExamples.C:119 IntervalExamples.C:120 IntervalExamples.C:121 IntervalExamples.C:122 IntervalExamples.C:123 IntervalExamples.C:124 IntervalExamples.C:125 IntervalExamples.C:126 IntervalExamples.C:127 IntervalExamples.C:128 IntervalExamples.C:129 IntervalExamples.C:130 IntervalExamples.C:131 IntervalExamples.C:132 IntervalExamples.C:133 IntervalExamples.C:134 IntervalExamples.C:135 IntervalExamples.C:136 IntervalExamples.C:137 IntervalExamples.C:138 IntervalExamples.C:139 IntervalExamples.C:140 IntervalExamples.C:141 IntervalExamples.C:142 IntervalExamples.C:143 IntervalExamples.C:144 IntervalExamples.C:145 IntervalExamples.C:146 IntervalExamples.C:147 IntervalExamples.C:148 IntervalExamples.C:149 IntervalExamples.C:150 IntervalExamples.C:151 IntervalExamples.C:152 IntervalExamples.C:153 IntervalExamples.C:154 IntervalExamples.C:155 IntervalExamples.C:156 IntervalExamples.C:157 IntervalExamples.C:158 IntervalExamples.C:159 IntervalExamples.C:160 IntervalExamples.C:161 IntervalExamples.C:162 IntervalExamples.C:163 IntervalExamples.C:164 IntervalExamples.C:165 IntervalExamples.C:166 IntervalExamples.C:167 IntervalExamples.C:168 IntervalExamples.C:169 IntervalExamples.C:170 IntervalExamples.C:171 IntervalExamples.C:172 IntervalExamples.C:173 IntervalExamples.C:174 IntervalExamples.C:175 IntervalExamples.C:176 IntervalExamples.C:177 IntervalExamples.C:178 IntervalExamples.C:179 IntervalExamples.C:180 IntervalExamples.C:181 IntervalExamples.C:182 IntervalExamples.C:183 IntervalExamples.C:184 IntervalExamples.C:185 IntervalExamples.C:186 IntervalExamples.C:187 IntervalExamples.C:188 IntervalExamples.C:189 IntervalExamples.C:190 IntervalExamples.C:191 IntervalExamples.C:192 IntervalExamples.C:193 IntervalExamples.C:194 IntervalExamples.C:195 IntervalExamples.C:196 IntervalExamples.C:197 IntervalExamples.C:198 IntervalExamples.C:199 IntervalExamples.C:200 IntervalExamples.C:201 IntervalExamples.C:202 IntervalExamples.C:203 IntervalExamples.C:204 IntervalExamples.C:205 IntervalExamples.C:206 IntervalExamples.C:207 IntervalExamples.C:208 IntervalExamples.C:209 IntervalExamples.C:210 IntervalExamples.C:211 IntervalExamples.C:212 IntervalExamples.C:213 IntervalExamples.C:214 IntervalExamples.C:215 IntervalExamples.C:216 IntervalExamples.C:217 IntervalExamples.C:218 |
|