TwoSidedFrequentistUpperLimitWithBands.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_roostats
/// \notebook -js
/// TwoSidedFrequentistUpperLimitWithBands
///
///
/// This is a standard demo that can be used with any ROOT file
/// prepared in the standard way.  You specify:
///  - name for input ROOT file
///  - name of workspace inside ROOT file that holds model and data
///  - name of ModelConfig that specifies details for calculator tools
///  - name of dataset
///
/// With default parameters the macro will attempt to run the
/// standard hist2workspace example and read the ROOT file
/// that it produces.
///
/// You may want to control:
/// ~~~{.cpp}
///   double confidenceLevel=0.95;
///   double additionalToysFac = 1.;
///   int nPointsToScan = 30;
///   int nToyMC = 500;
/// ~~~
///
/// This uses a modified version of the profile likelihood ratio as
/// a test statistic for upper limits (eg. test stat = 0 if muhat>mu).
///
/// Based on the observed data, one defines a set of parameter points
/// to be tested based on the value of the parameter of interest
/// and the conditional MLE (eg. profiled) values of the nuisance parameters.
///
/// At each parameter point, pseudo-experiments are generated using this
/// fixed reference model and then the test statistic is evaluated.
/// The auxiliary measurements (global observables) associated with the
/// constraint terms in nuisance parameters are also fluctuated in the
/// process of generating the pseudo-experiments in a frequentist manner
/// forming an 'unconditional ensemble'.  One could form a 'conditional'
/// ensemble in which these auxiliary measurements are fixed.  Note that the
/// nuisance parameters are not randomized, which is a Bayesian procedure.
/// Note, the nuisance parameters are floating in the fits.  For each point,
/// the threshold that defines the 95% acceptance region is found.  This
/// forms a "Confidence Belt".
///
/// After constructing the confidence belt, one can find the confidence
/// interval for any particular dataset by finding the intersection
/// of the observed test statistic and the confidence belt.  First
/// this is done on the observed data to get an observed 1-sided upper limt.
///
/// Finally, there expected limit and bands (from background-only) are
/// formed by generating background-only data and finding the upper limit.
/// The background-only is defined as such that the nuisance parameters are
/// fixed to their best fit value based on the data with the signal rate fixed to 0.
/// The bands are done by hand for now, will later be part of the RooStats tools.
///
/// On a technical note, this technique IS the generalization of Feldman-Cousins
/// with nuisance parameters.
///
/// Building the confidence belt can be computationally expensive.
/// Once it is built, one could save it to a file and use it in a separate step.
///
/// We can use PROOF to speed things along in parallel, however,
/// the test statistic has to be installed on the workers
/// so either turn off PROOF or include the modified test statistic
/// in your $ROOTSYS/roofit/roostats/inc directory,
/// add the additional line to the LinkDef.h file,
/// and recompile root.
///
/// Note, if you have a boundary on the parameter of interest (eg. cross-section)
/// the threshold on the two-sided test statistic starts off at moderate values and plateaus.
///
/// [#0] PROGRESS:Generation -- generated toys: 500 / 999
/// NeymanConstruction: Prog: 12/50 total MC = 39 this test stat = 0
///  SigXsecOverSM=0.69 alpha_syst1=0.136515 alpha_syst3=0.425415 beta_syst2=1.08496 [-1e+30, 0.011215]  in interval = 1
///
/// this tells you the values of the parameters being used to generate the pseudo-experiments
/// and the threshold in this case is 0.011215.  One would expect for 95% that the threshold
/// would be ~1.35 once the cross-section is far enough away from 0 that it is essentially
/// unaffected by the boundary.  As one reaches the last points in the scan, the
/// theshold starts to get artificially high.  This is because the range of the parameter in
/// the fit is the same as the range in the scan.  In the future, these should be independently
/// controlled, but they are not now.  As a result the ~50% of pseudo-experiments that have an
/// upward fluctuation end up with muhat = muMax.  Because of this, the upper range of the
/// parameter should be well above the expected upper limit... but not too high or one will
/// need a very large value of nPointsToScan to resolve the relevant region.  This can be
/// improved, but this is the first version of this script.
///
/// Important note: when the model includes external constraint terms, like a Gaussian
/// constraint to a nuisance parameter centered around some nominal value there is
/// a subtlety.  The asymptotic results are all based on the assumption that all the
/// measurements fluctuate... including the nominal values from auxiliary measurements.
/// If these do not fluctuate, this corresponds to an "conditional ensemble".  The
/// result is that the distribution of the test statistic can become very non-chi^2.
/// This results in thresholds that become very large.
///
/// \macro_image
/// \macro_output
/// \macro_code
///
/// \author Kyle Cranmer,Contributions from Aaron Armbruster, Haoshuang Ji, Haichen Wang and Daniel Whiteson

#include "TFile.h"
#include "TROOT.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TSystem.h"
#include <iostream>

#include "RooWorkspace.h"
#include "RooSimultaneous.h"
#include "RooAbsData.h"

#include "RooStats/ModelConfig.h"
#include "RooStats/FeldmanCousins.h"
#include "RooStats/ToyMCSampler.h"
#include "RooStats/PointSetInterval.h"
#include "RooStats/ConfidenceBelt.h"

#include "RooStats/RooStatsUtils.h"
#include "RooStats/ProfileLikelihoodTestStat.h"

using namespace RooFit;
using namespace RooStats;
using namespace std;

bool useProof = false;  // flag to control whether to use Proof
int nworkers = 0;   // number of workers (default use all available cores)

// -------------------------------------------------------

void TwoSidedFrequentistUpperLimitWithBands(const char* infile = "",
                                            const char* workspaceName = "combined",
                                            const char* modelConfigName = "ModelConfig",
                                            const char* dataName = "obsData") {


   double confidenceLevel=0.95;
   // degrade/improve number of pseudo-experiments used to define the confidence belt.
   // value of 1 corresponds to default number of toys in the tail, which is 50/(1-confidenceLevel)
   double additionalToysFac = 0.5;
   int nPointsToScan = 20; // number of steps in the parameter of interest
   int nToyMC = 200; // number of toys used to define the expected limit and band

   // -------------------------------------------------------
   // First part is just to access a user-defined file
   // or create the standard example file if it doesn't exist
   const char* filename = "";
      if (!strcmp(infile,"")) {
         filename = "results/example_combined_GaussExample_model.root";
         bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
         // if file does not exists generate with histfactory
         if (!fileExist) {
#ifdef _WIN32
            cout << "HistFactory file cannot be generated on Windows - exit" << endl;
            return;
#endif
            // Normally this would be run on the command line
            cout <<"will run standard hist2workspace example"<<endl;
            gROOT->ProcessLine(".! prepareHistFactory .");
            gROOT->ProcessLine(".! hist2workspace config/example.xml");
            cout <<"\n\n---------------------"<<endl;
            cout <<"Done creating example input"<<endl;
            cout <<"---------------------\n\n"<<endl;
         }

      }
      else
         filename = infile;

      // Try to open the file
      TFile *file = TFile::Open(filename);

      // if input file was specified byt not found, quit
      if(!file ){
         cout <<"StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
         return;
      }

   // -------------------------------------------------------
   // Now get the data and workspace

   // get the workspace out of the file
   RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);
   if(!w){
      cout <<"workspace not found" << endl;
      return;
   }

   // get the modelConfig out of the file
   ModelConfig* mc = (ModelConfig*) w->obj(modelConfigName);

   // get the modelConfig out of the file
   RooAbsData* data = w->data(dataName);

   // make sure ingredients are found
   if(!data || !mc){
      w->Print();
      cout << "data or ModelConfig was not found" <<endl;
      return;
   }

   cout << "Found data and ModelConfig:" <<endl;
   mc->Print();

   // -------------------------------------------------------
   // Now get the POI for convenience
   // you may want to adjust the range of your POI
   RooRealVar* firstPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
   /*  firstPOI->setMin(0);*/
   /*  firstPOI->setMax(10);*/

   // -------------------------------------------------------
   // create and use the FeldmanCousins tool
   // to find and plot the 95% confidence interval
   // on the parameter of interest as specified
   // in the model config
   // REMEMBER, we will change the test statistic
   // so this is NOT a Feldman-Cousins interval
   FeldmanCousins fc(*data,*mc);
   fc.SetConfidenceLevel(confidenceLevel);
   fc.AdditionalNToysFactor(additionalToysFac); // improve sampling that defines confidence belt
   //  fc.UseAdaptiveSampling(true); // speed it up a bit, but don't use for expected limits
   fc.SetNBins(nPointsToScan); // set how many points per parameter of interest to scan
   fc.CreateConfBelt(true); // save the information in the belt for plotting

   // -------------------------------------------------------
   // Feldman-Cousins is a unified limit by definition
   // but the tool takes care of a few things for us like which values
   // of the nuisance parameters should be used to generate toys.
   // so let's just change the test statistic and realize this is
   // no longer "Feldman-Cousins" but is a fully frequentist Neyman-Construction.
   //  fc.GetTestStatSampler()->SetTestStatistic(&onesided);
   // ((ToyMCSampler*) fc.GetTestStatSampler())->SetGenerateBinned(true);
   ToyMCSampler*  toymcsampler = (ToyMCSampler*) fc.GetTestStatSampler();
   ProfileLikelihoodTestStat* testStat = dynamic_cast<ProfileLikelihoodTestStat*>(toymcsampler->GetTestStatistic());

   // Since this tool needs to throw toy MC the PDF needs to be
   // extended or the tool needs to know how many entries in a dataset
   // per pseudo experiment.
   // In the 'number counting form' where the entries in the dataset
   // are counts, and not values of discriminating variables, the
   // datasets typically only have one entry and the PDF is not
   // extended.
   if(!mc->GetPdf()->canBeExtended()){
      if(data->numEntries()==1)
         fc.FluctuateNumDataEntries(false);
      else
         cout <<"Not sure what to do about this model" <<endl;
   }

   // We can use PROOF to speed things along in parallel
   // However, the test statistic has to be installed on the workers
   // so either turn off PROOF or include the modified test statistic
   // in your $ROOTSYS/roofit/roostats/inc directory,
   // add the additional line to the LinkDef.h file,
   // and recompile root.
   if (useProof) {
      ProofConfig pc(*w, nworkers, "",false);
      toymcsampler->SetProofConfig(&pc); // enable proof
   }

   if(mc->GetGlobalObservables()){
      cout << "will use global observables for unconditional ensemble"<<endl;
      mc->GetGlobalObservables()->Print();
      toymcsampler->SetGlobalObservables(*mc->GetGlobalObservables());
   }


   // Now get the interval
   PointSetInterval* interval = fc.GetInterval();
   ConfidenceBelt* belt = fc.GetConfidenceBelt();

   // print out the interval on the first Parameter of Interest
   cout << "\n95% interval on " <<firstPOI->GetName()<<" is : ["<<
      interval->LowerLimit(*firstPOI) << ", "<<
      interval->UpperLimit(*firstPOI) <<"] "<<endl;

   // get observed UL and value of test statistic evaluated there
   RooArgSet tmpPOI(*firstPOI);
   double observedUL = interval->UpperLimit(*firstPOI);
   firstPOI->setVal(observedUL);
   double obsTSatObsUL = fc.GetTestStatSampler()->EvaluateTestStatistic(*data,tmpPOI);


   // Ask the calculator which points were scanned
   RooDataSet* parameterScan = (RooDataSet*) fc.GetPointsToScan();
   RooArgSet* tmpPoint;

   // make a histogram of parameter vs. threshold
   TH1F* histOfThresholds = new TH1F("histOfThresholds","",
                                       parameterScan->numEntries(),
                                       firstPOI->getMin(),
                                       firstPOI->getMax());
   histOfThresholds->GetXaxis()->SetTitle(firstPOI->GetName());
   histOfThresholds->GetYaxis()->SetTitle("Threshold");

   // loop through the points that were tested and ask confidence belt
   // what the upper/lower thresholds were.
   // For FeldmanCousins, the lower cut off is always 0
   for(Int_t i=0; i<parameterScan->numEntries(); ++i){
      tmpPoint = (RooArgSet*) parameterScan->get(i)->clone("temp");
      //cout <<"get threshold"<<endl;
      double arMax = belt->GetAcceptanceRegionMax(*tmpPoint);
      double poiVal = tmpPoint->getRealValue(firstPOI->GetName()) ;
      histOfThresholds->Fill(poiVal,arMax);
   }
   TCanvas* c1 = new TCanvas();
   c1->Divide(2);
   c1->cd(1);
   histOfThresholds->SetMinimum(0);
   histOfThresholds->Draw();
   c1->cd(2);

   // -------------------------------------------------------
   // Now we generate the expected bands and power-constraint

   // First: find parameter point for mu=0, with conditional MLEs for nuisance parameters
   RooAbsReal* nll = mc->GetPdf()->createNLL(*data);
   RooAbsReal* profile = nll->createProfile(*mc->GetParametersOfInterest());
   firstPOI->setVal(0.);
   profile->getVal(); // this will do fit and set nuisance parameters to profiled values
   RooArgSet* poiAndNuisance = new RooArgSet();
   if(mc->GetNuisanceParameters())
      poiAndNuisance->add(*mc->GetNuisanceParameters());
   poiAndNuisance->add(*mc->GetParametersOfInterest());
   w->saveSnapshot("paramsToGenerateData",*poiAndNuisance);
   RooArgSet* paramsToGenerateData = (RooArgSet*) poiAndNuisance->snapshot();
   cout << "\nWill use these parameter points to generate pseudo data for bkg only" << endl;
   paramsToGenerateData->Print("v");


   RooArgSet unconditionalObs;
   unconditionalObs.add(*mc->GetObservables());
   unconditionalObs.add(*mc->GetGlobalObservables()); // comment this out for the original conditional ensemble

   double CLb=0;
   double CLbinclusive=0;

   // Now we generate background only and find distribution of upper limits
   TH1F* histOfUL = new TH1F("histOfUL","",100,0,firstPOI->getMax());
   histOfUL->GetXaxis()->SetTitle("Upper Limit (background only)");
   histOfUL->GetYaxis()->SetTitle("Entries");
   for(int imc=0; imc<nToyMC; ++imc){

      // set parameters back to values for generating pseudo data
      //    cout << "\n get current nuis, set vals, print again" << endl;
      w->loadSnapshot("paramsToGenerateData");
      //    poiAndNuisance->Print("v");

      RooDataSet* toyData = 0;
      // now generate a toy dataset for the main measurement
      if(!mc->GetPdf()->canBeExtended()){
         if(data->numEntries()==1)
            toyData = mc->GetPdf()->generate(*mc->GetObservables(),1);
         else
            cout <<"Not sure what to do about this model" <<endl;
      } else{
         //      cout << "generating extended dataset"<<endl;
         toyData = mc->GetPdf()->generate(*mc->GetObservables(),Extended());
      }

      // generate global observables
      // need to be careful for simpdf.
      // In ROOT 5.28 there is a problem with generating global observables
      // with a simultaneous PDF.  In 5.29 there is a solution with
      // RooSimultaneous::generateSimGlobal, but this may change to
      // the standard generate interface in 5.30.

      RooSimultaneous* simPdf = dynamic_cast<RooSimultaneous*>(mc->GetPdf());
      if(!simPdf){
         RooDataSet *one = mc->GetPdf()->generate(*mc->GetGlobalObservables(), 1);
         const RooArgSet *values = one->get();
         RooArgSet *allVars = mc->GetPdf()->getVariables();
         *allVars = *values;
         delete allVars;
         delete one;
      } else {
         RooDataSet* one = simPdf->generateSimGlobal(*mc->GetGlobalObservables(),1);
         const RooArgSet *values = one->get();
         RooArgSet *allVars = mc->GetPdf()->getVariables();
         *allVars = *values;
         delete allVars;
         delete one;

      }


      // get test stat at observed UL in observed data
      firstPOI->setVal(observedUL);
      double toyTSatObsUL = fc.GetTestStatSampler()->EvaluateTestStatistic(*toyData,tmpPOI);
      //    toyData->get()->Print("v");
      //    cout <<"obsTSatObsUL " <<obsTSatObsUL << "toyTS " << toyTSatObsUL << endl;
      if(obsTSatObsUL < toyTSatObsUL) // not sure about <= part yet
         CLb+= (1.)/nToyMC;
      if(obsTSatObsUL <= toyTSatObsUL) // not sure about <= part yet
         CLbinclusive+= (1.)/nToyMC;


      // loop over points in belt to find upper limit for this toy data
      double thisUL = 0;
      for(Int_t i=0; i<parameterScan->numEntries(); ++i){
         tmpPoint = (RooArgSet*) parameterScan->get(i)->clone("temp");
         double arMax = belt->GetAcceptanceRegionMax(*tmpPoint);
         firstPOI->setVal( tmpPoint->getRealValue(firstPOI->GetName()) );
         //   double thisTS = profile->getVal();
         double thisTS = fc.GetTestStatSampler()->EvaluateTestStatistic(*toyData,tmpPOI);

         //   cout << "poi = " << firstPOI->getVal()
         // << " max is " << arMax << " this profile = " << thisTS << endl;
         //      cout << "thisTS = " << thisTS<<endl;
         if(thisTS<=arMax){
            thisUL = firstPOI->getVal();
         } else{
            break;
         }
      }


      histOfUL->Fill(thisUL);

      // for few events, data is often the same, and UL is often the same
      //    cout << "thisUL = " << thisUL<<endl;

      delete toyData;
   }
   histOfUL->Draw();
   c1->SaveAs("two-sided_upper_limit_output.pdf");

   // if you want to see a plot of the sampling distribution for a particular scan point:
   /*
   SamplingDistPlot sampPlot;
   int indexInScan = 0;
   tmpPoint = (RooArgSet*) parameterScan->get(indexInScan)->clone("temp");
   firstPOI->setVal( tmpPoint->getRealValue(firstPOI->GetName()) );
   toymcsampler->SetParametersForTestStat(tmpPOI);
   SamplingDistribution* samp = toymcsampler->GetSamplingDistribution(*tmpPoint);
   sampPlot.AddSamplingDistribution(samp);
   sampPlot.Draw();
      */

   // Now find bands and power constraint
   Double_t* bins = histOfUL->GetIntegral();
   TH1F* cumulative = (TH1F*) histOfUL->Clone("cumulative");
   cumulative->SetContent(bins);
   double band2sigDown=0, band1sigDown=0, bandMedian=0, band1sigUp=0,band2sigUp=0;
   for(int i=1; i<=cumulative->GetNbinsX(); ++i){
      if(bins[i]<RooStats::SignificanceToPValue(2))
         band2sigDown=cumulative->GetBinCenter(i);
      if(bins[i]<RooStats::SignificanceToPValue(1))
         band1sigDown=cumulative->GetBinCenter(i);
      if(bins[i]<0.5)
         bandMedian=cumulative->GetBinCenter(i);
      if(bins[i]<RooStats::SignificanceToPValue(-1))
         band1sigUp=cumulative->GetBinCenter(i);
      if(bins[i]<RooStats::SignificanceToPValue(-2))
         band2sigUp=cumulative->GetBinCenter(i);
   }
   cout << "-2 sigma  band " << band2sigDown << endl;
   cout << "-1 sigma  band " << band1sigDown << " [Power Constraint)]" << endl;
   cout << "median of band " << bandMedian << endl;
   cout << "+1 sigma  band " << band1sigUp << endl;
   cout << "+2 sigma  band " << band2sigUp << endl;

   // print out the interval on the first Parameter of Interest
   cout << "\nobserved 95% upper-limit "<< interval->UpperLimit(*firstPOI) <<endl;
   cout << "CLb strict [P(toy>obs|0)] for observed 95% upper-limit "<< CLb <<endl;
   cout << "CLb inclusive [P(toy>=obs|0)] for observed 95% upper-limit "<< CLbinclusive <<endl;

   delete profile;
   delete nll;

}