HypoTestResult.cxx

Go to the documentation of this file.

// @(#)root/roostats:$Id$
// Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke, Sven Kreiss
/*************************************************************************
 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/
 
/*****************************************************************************
 * Project: RooStats
 * Package: RooFit/RooStats
 * @(#)root/roofit/roostats:$Id$
 * Authors:
 *   Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke, Sven Kreiss
 *
 *****************************************************************************/
 
 
/** \class RooStats::HypoTestResult
    \ingroup Roostats
 
HypoTestResult is a base class for results from hypothesis tests.
Any tool inheriting from HypoTestCalculator can return a HypoTestResult.
As such, it stores a p-value for the null-hypothesis (eg. background-only)
and an alternate hypothesis (eg. signal+background).
The p-values can also be transformed into confidence levels
(\f$CL_{b}\f$, \f$CL_{s+b}\f$) in a trivial way.
The ratio of the \f$CL_{s+b}\f$ to \f$CL_{b}\f$ is often called
\f$CL_{s}\f$, and is considered useful, though it is not a probability.
Finally, the p-value of the null can be transformed into a number of
equivalent Gaussian sigma using the Significance method.
 
The p-value of the null for a given test statistic is rigorously defined and
this is the starting point for the following conventions.
 
### Conventions used in this class
 
The p-value for the null and alternate are on the **same side** of the
observed value of the test statistic. This is the more standard
convention and avoids confusion when doing inverted tests.
 
For exclusion, we also want the formula \f$CL_{s} = CL_{s+b} / CL_{b}\f$
to hold which therefore defines our conventions for \f$CL_{s+b}\f$ and
\f$CL_{b}\f$. \f$CL_{s}\f$ was specifically invented for exclusion
and therefore all quantities need be related through the assignments
as they are for exclusion: \f$CL_{s+b} = p_{s+b}\f$; \f$CL_{b} = p_{b}\f$. This
is derived by considering the scenarios of a powerful and not powerful
inverted test, where for the not so powerful test, \f$CL_{s}\f$ must be
close to one.
 
For results of Hypothesis tests,
\f$CL_{s}\f$ has no similar direct interpretation as for exclusion and can
be larger than one.
 
*/
 
#include "RooStats/HypoTestResult.h"
#include "RooStats/SamplingDistribution.h"
#include "RooAbsReal.h"
 
#include "RooStats/RooStatsUtils.h"
 
#include <TMath.h>
 
#include <limits>
#define NaN numeric_limits<float>::quiet_NaN()
#define IsNaN(a) TMath::IsNaN(a)
 
 
using namespace RooStats;
using std::numeric_limits, std::endl;
 
////////////////////////////////////////////////////////////////////////////////
/// Default constructor
 
HypoTestResult::HypoTestResult(const char* name) :
   TNamed(name,name),
   fNullPValue(NaN), fAlternatePValue(NaN),
   fNullPValueError(0), fAlternatePValueError(0),
   fTestStatisticData(NaN),
   fAllTestStatisticsData(nullptr),
   fNullDistr(nullptr), fAltDistr(nullptr),
   fNullDetailedOutput(nullptr), fAltDetailedOutput(nullptr),
   fPValueIsRightTail(true),
   fBackgroundIsAlt(false)
{
}
 
////////////////////////////////////////////////////////////////////////////////
/// Alternate constructor
 
HypoTestResult::HypoTestResult(const char* name, double nullp, double altp) :
   TNamed(name,name),
   fNullPValue(nullp), fAlternatePValue(altp),
   fNullPValueError(0), fAlternatePValueError(0),
   fTestStatisticData(NaN),
   fAllTestStatisticsData(nullptr),
   fNullDistr(nullptr), fAltDistr(nullptr),
   fNullDetailedOutput(nullptr), fAltDetailedOutput(nullptr),
   fPValueIsRightTail(true),
   fBackgroundIsAlt(false)
{
}
 
////////////////////////////////////////////////////////////////////////////////
/// copy constructor
 
HypoTestResult::HypoTestResult(const HypoTestResult& other) :
   TNamed(other),
   fNullPValue(other.fNullPValue), fAlternatePValue(other.fAlternatePValue),
   fNullPValueError(other.fNullPValueError), fAlternatePValueError(other.fAlternatePValueError),
   fTestStatisticData(other.fTestStatisticData),
   fAllTestStatisticsData(nullptr),
   fNullDistr(nullptr), fAltDistr(nullptr),
   fNullDetailedOutput(nullptr), fAltDetailedOutput(nullptr),
   fPValueIsRightTail( other.GetPValueIsRightTail() ),
   fBackgroundIsAlt( other.GetBackGroundIsAlt() )
{
   this->Append( &other );
}
 
////////////////////////////////////////////////////////////////////////////////
/// Destructor
 
HypoTestResult::~HypoTestResult()
{
 
}
 
////////////////////////////////////////////////////////////////////////////////
/// assignment operator
 
HypoTestResult & HypoTestResult::operator=(const HypoTestResult& other) {
   if (this == &other) return *this;
   SetName(other.GetName());
   SetTitle(other.GetTitle());
   fNullPValue = other.fNullPValue;
   fAlternatePValue = other.fAlternatePValue;
   fNullPValueError = other.fNullPValueError;
   fAlternatePValueError = other.fAlternatePValueError;
   fTestStatisticData = other.fTestStatisticData;
 
 
   fAllTestStatisticsData = nullptr;
   fNullDistr = nullptr;
   fAltDistr = nullptr;
   fNullDetailedOutput = nullptr;
   fAltDetailedOutput = nullptr;
 
   fFitInfo = nullptr;
 
   fPValueIsRightTail =  other.GetPValueIsRightTail();
   fBackgroundIsAlt = other.GetBackGroundIsAlt();
 
   this->Append( &other );
 
   return *this;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Add additional toy-MC experiments to the current results.
/// Use the data test statistics of the added object if it is not already
/// set (otherwise, ignore the new one).
 
void HypoTestResult::Append(const HypoTestResult* other) {
   if (fNullDistr) {
      fNullDistr->Add(other->GetNullDistribution());
   } else if (other->GetNullDistribution()) {
      fNullDistr = std::make_unique<SamplingDistribution>(*other->GetNullDistribution());
   }
 
   if (fAltDistr) {
      fAltDistr->Add(other->GetAltDistribution());
   } else if (other->GetAltDistribution()) {
      fAltDistr = std::make_unique<SamplingDistribution>(*other->GetAltDistribution());
   }
 
   if( fNullDetailedOutput ) {
      if( other->GetNullDetailedOutput() ) fNullDetailedOutput->append( *other->GetNullDetailedOutput() );
   }else{
      if( other->GetNullDetailedOutput() ) fNullDetailedOutput = std::make_unique<RooDataSet>( *other->GetNullDetailedOutput() );
   }
 
   if( fAltDetailedOutput ) {
      if( other->GetAltDetailedOutput() ) fAltDetailedOutput->append( *other->GetAltDetailedOutput() );
   }else{
      if( other->GetAltDetailedOutput() ) fAltDetailedOutput = std::make_unique<RooDataSet>( *other->GetAltDetailedOutput() );
   }
 
   if( fFitInfo ) {
      if( other->GetFitInfo() ) fFitInfo->append( *other->GetFitInfo() );
   }else{
      if( other->GetFitInfo() ) fFitInfo = std::make_unique<RooDataSet>( *other->GetFitInfo());
   }
 
   // if no data is present use the other HypoTestResult's data
   if(IsNaN(fTestStatisticData)) fTestStatisticData = other->GetTestStatisticData();
 
   UpdatePValue(fNullDistr.get(), fNullPValue, fNullPValueError, true);
   UpdatePValue(fAltDistr.get(), fAlternatePValue, fAlternatePValueError, false);
}
 
////////////////////////////////////////////////////////////////////////////////
 
void HypoTestResult::SetAltDistribution(SamplingDistribution *alt) {
   fAltDistr.reset(alt);
   UpdatePValue(fAltDistr.get(), fAlternatePValue, fAlternatePValueError, false);
}
 
////////////////////////////////////////////////////////////////////////////////
 
void HypoTestResult::SetNullDistribution(SamplingDistribution *null) {
   fNullDistr.reset(null);
   UpdatePValue(fNullDistr.get(), fNullPValue, fNullPValueError, true);
}
 
////////////////////////////////////////////////////////////////////////////////
 
void HypoTestResult::SetTestStatisticData(const double tsd) {
   fTestStatisticData = tsd;
 
   UpdatePValue(fNullDistr.get(), fNullPValue, fNullPValueError, true);
   UpdatePValue(fAltDistr.get(), fAlternatePValue, fAlternatePValueError, false);
}
 
////////////////////////////////////////////////////////////////////////////////
 
void HypoTestResult::SetAllTestStatisticsData(const RooArgList* tsd) {
   if (tsd) fAllTestStatisticsData.reset(static_cast<const RooArgList*>(tsd->snapshot()));
 
   if( fAllTestStatisticsData  &&  !fAllTestStatisticsData->empty() ) {
      RooRealVar* firstTS = static_cast<RooRealVar*>(fAllTestStatisticsData->at(0));
      if( firstTS ) SetTestStatisticData( firstTS->getVal() );
   }
}
 
////////////////////////////////////////////////////////////////////////////////
 
void HypoTestResult::SetPValueIsRightTail(bool pr) {
   fPValueIsRightTail = pr;
 
   UpdatePValue(fNullDistr.get(), fNullPValue, fNullPValueError, true);
   UpdatePValue(fAltDistr.get(), fAlternatePValue, fAlternatePValueError, false);
}
 
////////////////////////////////////////////////////////////////////////////////
 
bool HypoTestResult::HasTestStatisticData(void) const {
   return !IsNaN(fTestStatisticData);
}
 
////////////////////////////////////////////////////////////////////////////////
 
double HypoTestResult::NullPValueError() const {
   // compute error on Null pvalue
   return fNullPValueError;
}
 
////////////////////////////////////////////////////////////////////////////////
/// compute \f$CL_{b}\f$ error
/// \f$CL_{b}\f$ = 1 - NullPValue()
/// must use opposite condition that routine above
 
double HypoTestResult::CLbError() const {
   return fBackgroundIsAlt ? fAlternatePValueError : fNullPValueError;
}
 
////////////////////////////////////////////////////////////////////////////////
 
double HypoTestResult::CLsplusbError() const {
   return fBackgroundIsAlt ? fNullPValueError : fAlternatePValueError;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Taylor expansion series approximation for standard deviation (error propagation)
 
double HypoTestResult::SignificanceError() const {
   return NullPValueError() / ROOT::Math::normal_pdf(Significance());
}
 
////////////////////////////////////////////////////////////////////////////////
/// Returns an estimate of the error on \f$CL_{s}\f$ through combination of the
/// errors on \f$CL_{b}\f$ and \f$CL_{s+b}\f$:
/// \f[
/// \sigma_{CL_s} = CL_s
/// \sqrt{\left( \frac{\sigma_{CL_{s+b}}}{CL_{s+b}} \right)^2 + \left( \frac{\sigma_{CL_{b}}}{CL_{b}} \right)^2}
/// \f]
 
double HypoTestResult::CLsError() const {
   if(!fAltDistr || !fNullDistr) return 0.0;
 
   // unsigned const int n_b = fNullDistr->GetSamplingDistribution().size();
   // unsigned const int n_sb = fAltDistr->GetSamplingDistribution().size();
 
   // if CLb() == 0 CLs = -1 so return a -1 error
   if (CLb() == 0 ) return -1;
 
   double cl_b_err2 = pow(CLbError(),2);
   double cl_sb_err2 = pow(CLsplusbError(),2);
 
   return std::sqrt(cl_sb_err2 + cl_b_err2 * pow(CLs(),2))/CLb();
}
 
////////////////////////////////////////////////////////////////////////////////
/// updates the pvalue if sufficient data is available
 
void HypoTestResult::UpdatePValue(const SamplingDistribution* distr, double &pvalue, double &perror, bool /*isNull*/) {
   if(IsNaN(fTestStatisticData)) return;
   if(!distr) return;
 
   /* Got to be careful for discrete distributions:
    * To get the right behaviour for limits, the p-value must
    * include the value of fTestStatistic both for Alt and Null cases
    */
   if(fPValueIsRightTail) {
      pvalue = distr->IntegralAndError(perror, fTestStatisticData, RooNumber::infinity(), true,
                                       true , true );   // always closed interval [ fTestStatistic, inf ]
 
   }else{
      pvalue = distr->IntegralAndError(perror, -RooNumber::infinity(), fTestStatisticData, true,
                                       true,  true  ); // // always closed  [ -inf, fTestStatistic ]
   }
}
 
////////////////////////////////////////////////////////////////////////////////
/// Print out some information about the results
/// Note: use Alt/Null labels for the hypotheses here as the Null
/// might be the s+b hypothesis.
 
void HypoTestResult::Print(Option_t * ) const
{
   bool fromToys = (fAltDistr || fNullDistr);
 
   std::cout << std::endl << "Results " << GetName() << ": " << std::endl;
   std::cout << " - Null p-value = " << NullPValue();
   if (fromToys) std::cout << " +/- " << NullPValueError();
   std::cout << std::endl;
   std::cout << " - Significance = " << Significance();
   if (fromToys) std::cout << " +/- " << SignificanceError() << " sigma";
   std::cout << std::endl;
   if(fAltDistr)
      std::cout << " - Number of Alt toys: " << fAltDistr->GetSize() << std::endl;
   if(fNullDistr)
      std::cout << " - Number of Null toys: " << fNullDistr->GetSize() << std::endl;
 
   if (HasTestStatisticData() ) std::cout << " - Test statistic evaluated on data: " << fTestStatisticData << std::endl;
   std::cout << " - CL_b: " << CLb();
   if (fromToys) std::cout << " +/- " << CLbError();
   std::cout << std::endl;
   std::cout << " - CL_s+b: " << CLsplusb();
   if (fromToys) std::cout << " +/- " << CLsplusbError();
   std::cout << std::endl;
   std::cout << " - CL_s: " << CLs();
   if (fromToys) std::cout << " +/- " << CLsError();
   std::cout << std::endl;
 
   return;
}