Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 31 additions & 0 deletions roofit/roostats/inc/RooStats/AsymptoticCalculator.h
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,11 @@
#include "RooArgSet.h"
#include "Rtypes.h"

#include <memory>

class RooArgList;
class RooCategory;
class RooFitResult;
class RooRealVar;
class RooPoisson;
class RooProdPdf;
Expand Down Expand Up @@ -93,6 +96,29 @@ namespace RooStats {
/// return best fit value for all parameters
const RooArgSet & GetBestFitParams() const { return fBestFitPoi; }

/// Result of the unconditional fit to the observed data, performed by
/// Initialize() (updated if GetHypoTest() finds a better minimum).
/// Returns nullptr if the fit was skipped or has not been run yet.
/// The calculator keeps ownership of the returned object.
const RooFitResult *GetFitResultUncondObs() const { return fFitResultUncondObs.get(); }
/// Result of the conditional fit to the observed data with the POI fixed
/// to the tested value, from the last call to GetHypoTest().
/// Returns nullptr if the fit was skipped or has not been run yet.
/// The calculator keeps ownership of the returned object.
const RooFitResult *GetFitResultCondObs() const { return fFitResultCondObs.get(); }
/// Result of the fit to the Asimov data set with the POI fixed to the
/// value of the alternate-model snapshot, performed by Initialize(). Since
/// the Asimov data set is generated at that POI value, this corresponds to
/// the unconditional minimum (updated if GetHypoTest() finds a better
/// minimum). Returns nullptr if the fit was skipped or has not been run
/// yet. The calculator keeps ownership of the returned object.
const RooFitResult *GetFitResultUncondAsimov() const { return fFitResultUncondAsimov.get(); }
/// Result of the conditional fit to the Asimov data set with the POI
/// fixed to the tested value, from the last call to GetHypoTest().
/// Returns nullptr if the fit was skipped or has not been run yet.
/// The calculator keeps ownership of the returned object.
const RooFitResult *GetFitResultCondAsimov() const { return fFitResultCondAsimov.get(); }

static void SetPrintLevel(int level);

private:
Expand All @@ -109,6 +135,11 @@ namespace RooStats {
mutable RooArgSet fBestFitPoi; ///< snapshot of best fitted POI values
mutable RooArgSet fBestFitParams; ///< snapshot of all best fitted Parameter values

mutable std::unique_ptr<RooFitResult> fFitResultUncondObs; ///<! result of unconditional fit to observed data
mutable std::unique_ptr<RooFitResult> fFitResultCondObs; ///<! result of conditional fit to observed data
mutable std::unique_ptr<RooFitResult> fFitResultUncondAsimov; ///<! result of fit to Asimov data at the alt POI
mutable std::unique_ptr<RooFitResult> fFitResultCondAsimov; ///<! result of conditional fit to Asimov data

ClassDefOverride(AsymptoticCalculator,0)
};
}
Expand Down
83 changes: 66 additions & 17 deletions roofit/roostats/src/AsymptoticCalculator.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,18 @@ If more than one POI exists, only the first one is used.
The calculator can generate Asimov datasets from two kinds of PDFs:
- "Counting" distributions: RooPoisson, RooGaussian, or products of RooPoissons.
- Extended, *i.e.* number of events can be read off from extended likelihood term.
*/

The fits performed by the calculator can be steered with the global default
minimizer options, *e.g.* via ROOT::Math::MinimizerOptions::SetDefaultStrategy()
and ROOT::Math::MinimizerOptions::SetDefaultTolerance() (the tolerance is
clamped to a minimum value of 1). The RooFitResult objects of the fits are
retrievable after calling GetHypoTest() via GetFitResultUncondObs(),
GetFitResultCondObs(), GetFitResultUncondAsimov() and GetFitResultCondAsimov(),
so quantities like the minimizer status or the EDM at the minimum can be
inspected, for example to cross-check a fit that did not converge. Note that
when the calculator is driven by the HypoTestInverter, the stored conditional
fit results correspond to the last scanned point.
*/

#include "RooStats/AsymptoticCalculator.h"
#include "RooStats/ModelConfig.h"
Expand Down Expand Up @@ -86,7 +96,8 @@ int &fgPrintLevel()
}

// Forward declaration.
double EvaluateNLL(RooStats::ModelConfig const &modelConfig, RooAbsData &data, const RooArgSet *poiSet = nullptr);
double EvaluateNLL(RooStats::ModelConfig const &modelConfig, RooAbsData &data, const RooArgSet *poiSet = nullptr,
std::unique_ptr<RooFitResult> *fitResult = nullptr);

} // namespace

Expand Down Expand Up @@ -198,11 +209,19 @@ bool AsymptoticCalculator::Initialize() const {
fBestFitPoi.removeAll();
fBestFitParams.removeAll();
fAsimovGlobObs.removeAll();
fFitResultUncondObs.reset();
fFitResultCondObs.reset();
fFitResultUncondAsimov.reset();
fFitResultCondAsimov.reset();

// evaluate the unconditional nll for the full model on the observed data
if (verbose >= 0)
oocoutP(nullptr,Eval) << "AsymptoticCalculator::Initialize - Find best unconditional NLL on observed data" << std::endl;
fNLLObs = EvaluateNLL(*GetNullModel(), data);
fNLLObs = EvaluateNLL(*GetNullModel(), data, nullptr, &fFitResultUncondObs);
if (fFitResultUncondObs) {
fFitResultUncondObs->SetName("fitResultUncondObs");
fFitResultUncondObs->SetTitle("Unconditional fit to observed data");
}
// fill also snapshot of best poi
poi->snapshot(fBestFitPoi);
RooRealVar * muBest = dynamic_cast<RooRealVar*>(fBestFitPoi.first());
Expand Down Expand Up @@ -285,7 +304,11 @@ bool AsymptoticCalculator::Initialize() const {
<< muAlt->GetName() << " ) = " << muAlt->getVal() << std::endl;
}

fNLLAsimov = EvaluateNLL(*GetNullModel(), *fAsimovData, &poiAlt );
fNLLAsimov = EvaluateNLL(*GetNullModel(), *fAsimovData, &poiAlt, &fFitResultUncondAsimov);
if (fFitResultUncondAsimov) {
fFitResultUncondAsimov->SetName("fitResultUncondAsimov");
fFitResultUncondAsimov->SetTitle("Fit to Asimov data with POI fixed to the alt-model snapshot");
}
// for unconditional fit
//fNLLAsimov = EvaluateNLL( *nullPdf, *fAsimovData);
//poi->Print("v");
Expand All @@ -302,10 +325,14 @@ bool AsymptoticCalculator::Initialize() const {

namespace {

double EvaluateNLL(RooStats::ModelConfig const& modelConfig, RooAbsData& data, const RooArgSet *poiSet)
double EvaluateNLL(RooStats::ModelConfig const &modelConfig, RooAbsData &data, const RooArgSet *poiSet,
std::unique_ptr<RooFitResult> *fitResult)
{
int verbose = fgPrintLevel();

if (fitResult)
fitResult->reset();

RooAbsPdf &pdf = *modelConfig.GetPdf();

RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
Expand Down Expand Up @@ -417,13 +444,12 @@ double EvaluateNLL(RooStats::ModelConfig const& modelConfig, RooAbsData& data, c
}
}

std::unique_ptr<RooFitResult> result;
// save the fit result also in case of failure, so that the status of a
// non-converged fit can be inspected by the user
std::unique_ptr<RooFitResult> result{minim.save()};

// ignore errors in Hesse or in Improve and also when matrix was made pos def (status returned = 1)
if (status >= 0) {
result = std::unique_ptr<RooFitResult>{minim.save()};
}
if (result){
if (status >= 0 && result) {
if (RooStats::NLLOffsetMode() != "initial") {
val = result->minNll();
} else {
Expand All @@ -433,12 +459,14 @@ double EvaluateNLL(RooStats::ModelConfig const& modelConfig, RooAbsData& data, c
if (!previous) RooAbsReal::setHideOffset(false) ;
}

}
else {
} else {
oocoutE(nullptr,Fitting) << "FIT FAILED !- return a NaN NLL " << std::endl;
val = TMath::QuietNaN();
}

if (fitResult)
*fitResult = std::move(result);

minim.optimizeConst(false);
}

Expand Down Expand Up @@ -530,7 +558,12 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const {
}

// evaluate the conditional NLL on the observed data for the snapshot value
double condNLL = EvaluateNLL(*GetNullModel(), const_cast<RooAbsData&>(*GetData()), &poiTest);
double condNLL = EvaluateNLL(*GetNullModel(), const_cast<RooAbsData &>(*GetData()), &poiTest, &fFitResultCondObs);
if (fFitResultCondObs) {
fFitResultCondObs->SetName("fitResultCondObs");
fFitResultCondObs->SetTitle(
TString::Format("Conditional fit to observed data for %s = %g", muTest->GetName(), muTest->getVal()));
}

double qmu = 2.*(condNLL - fNLLObs);

Expand All @@ -552,14 +585,20 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const {
<< "AsymptoticCalculator: unconditional fit failed before - retry to do it now " << std::endl;
}

double nll = EvaluateNLL(*GetNullModel(), const_cast<RooAbsData&>(*GetData()));
std::unique_ptr<RooFitResult> refitResult;
double nll = EvaluateNLL(*GetNullModel(), const_cast<RooAbsData &>(*GetData()), nullptr, &refitResult);

if (nll < fNLLObs || (TMath::IsNaN(fNLLObs) && !TMath::IsNaN(nll) ) ) {
oocoutW(nullptr,Minimization) << "AsymptoticCalculator: Found a better unconditional minimum "
<< " old NLL = " << fNLLObs << " old muHat " << muHat->getVal() << std::endl;

// update values
fNLLObs = nll;
if (refitResult) {
fFitResultUncondObs = std::move(refitResult);
fFitResultUncondObs->SetName("fitResultUncondObs");
fFitResultUncondObs->SetTitle("Unconditional fit to observed data");
}
const RooArgSet * poi = GetNullModel()->GetParametersOfInterest();
assert(poi);
fBestFitPoi.removeAll();
Expand Down Expand Up @@ -614,8 +653,12 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const {

if (verbose > 0) oocoutP(nullptr,Eval) << "AsymptoticCalculator::GetHypoTest -- Find best conditional NLL on ASIMOV data set .... " << std::endl;

double condNLL_A = EvaluateNLL(*GetNullModel(), *fAsimovData, &poiTest);

double condNLL_A = EvaluateNLL(*GetNullModel(), *fAsimovData, &poiTest, &fFitResultCondAsimov);
if (fFitResultCondAsimov) {
fFitResultCondAsimov->SetName("fitResultCondAsimov");
fFitResultCondAsimov->SetTitle(
TString::Format("Conditional fit to Asimov data for %s = %g", muTest->GetName(), muTest->getVal()));
}

double qmu_A = 2.*(condNLL_A - fNLLAsimov );

Expand All @@ -634,14 +677,20 @@ HypoTestResult* AsymptoticCalculator::GetHypoTest() const {
<< std::endl;
}

double nll = EvaluateNLL(*GetNullModel(), *fAsimovData);
std::unique_ptr<RooFitResult> refitResult;
double nll = EvaluateNLL(*GetNullModel(), *fAsimovData, nullptr, &refitResult);

if (nll < fNLLAsimov || (TMath::IsNaN(fNLLAsimov) && !TMath::IsNaN(nll) )) {
oocoutW(nullptr,Minimization) << "AsymptoticCalculator: Found a better unconditional minimum for Asimov data set"
<< " old NLL = " << fNLLAsimov << std::endl;

// update values
fNLLAsimov = nll;
if (refitResult) {
fFitResultUncondAsimov = std::move(refitResult);
fFitResultUncondAsimov->SetName("fitResultUncondAsimov");
fFitResultUncondAsimov->SetTitle("Unconditional fit to Asimov data");
}

oocoutW(nullptr,Minimization) << "AsymptoticCalculator: New minimum found for "
<< " NLL = " << fNLLAsimov << std::endl;
Expand Down
2 changes: 1 addition & 1 deletion roofit/roostats/test/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
ROOT_ADD_GTEST(testAsymptoticCalculator testAsymptoticCalculator.cxx LIBRARIES RooStats)
ROOT_ADD_GTEST(testAsymptoticCalculator testAsymptoticCalculator.cxx LIBRARIES RooStats ROOT::TestSupport)
ROOT_ADD_GTEST(testBayesianCalculator testBayesianCalculator.cxx LIBRARIES RooStats)
ROOT_ADD_GTEST(testHypoTestInvResult testHypoTestInvResult.cxx
LIBRARIES RooStats
Expand Down
88 changes: 88 additions & 0 deletions roofit/roostats/test/testAsymptoticCalculator.cxx
Original file line number Diff line number Diff line change
@@ -1,7 +1,15 @@
// Author: Jonas Rembser, CERN 01/2025

#include "RooDataSet.h"
#include "RooFitResult.h"
#include "RooMultiVarGaussian.h"
#include "RooRealVar.h"
#include "RooWorkspace.h"
#include "RooStats/AsymptoticCalculator.h"
#include "RooStats/HypoTestResult.h"
#include "RooStats/ModelConfig.h"

#include "ROOT/TestSupport.hxx"

#include "gtest/gtest.h"

Expand Down Expand Up @@ -37,3 +45,83 @@ TEST(AsymptoticCalculator, CountingAsimovDataSetFromMultiVarGaussian)
EXPECT_EQ(dataX.getVal(), mu.getVal());
}
}

// Check that the fit results of the fits performed in Initialize() and
// GetHypoTest() are stored and can be accessed by the user (JIRA ROOT-10066).
TEST(AsymptoticCalculator, StoredFitResults)
{
using namespace RooStats;

// On/off Poisson counting model: signal region with s + b expected events,
// control region constraining the background via tau * b.
RooWorkspace ws;
ws.factory("Poisson::px(x[150, 0, 500], sum::splusb(s[0, 0, 100], b[100, 0, 300]))");
ws.factory("Poisson::py(y[100, 0, 500], prod::taub(tau[1.0], b))");
ws.factory("PROD::model(px, py)");

RooRealVar &s = *ws.var("s");
RooArgSet obs{*ws.var("x"), *ws.var("y")};

RooDataSet data{"data", "data", obs};
data.add(obs);

ModelConfig sbModel{"sbModel", &ws};
sbModel.SetPdf(*ws.pdf("model"));
sbModel.SetObservables(obs);
sbModel.SetParametersOfInterest(RooArgSet{s});
sbModel.SetNuisanceParameters(RooArgSet{*ws.var("b")});
s.setVal(50.0);
sbModel.SetSnapshot(RooArgSet{s});

std::unique_ptr<ModelConfig> bModel{static_cast<ModelConfig *>(sbModel.Clone("bModel"))};
s.setVal(0.0);
bModel->SetSnapshot(RooArgSet{s});

// Some of the fits start at parameter values that already correspond to the
// minimum, in which case Minuit2 emits a harmless line-search warning.
ROOT::TestSupport::CheckDiagsRAII checkDiag;
checkDiag.optionalDiag(kWarning, "Minuit2", "VariableMetricBuilder No improvement in line search", false);

AsymptoticCalculator calc{data, *bModel, sbModel};
calc.SetOneSided(true);

// Before running the hypothesis test, no fit results are available.
EXPECT_EQ(calc.GetFitResultCondObs(), nullptr);
EXPECT_EQ(calc.GetFitResultCondAsimov(), nullptr);

std::unique_ptr<HypoTestResult> result{calc.GetHypoTest()};
ASSERT_NE(result, nullptr);

const RooFitResult *uncondObs = calc.GetFitResultUncondObs();
const RooFitResult *condObs = calc.GetFitResultCondObs();
const RooFitResult *uncondAsimov = calc.GetFitResultUncondAsimov();
const RooFitResult *condAsimov = calc.GetFitResultCondAsimov();

ASSERT_NE(uncondObs, nullptr);
ASSERT_NE(condObs, nullptr);
ASSERT_NE(uncondAsimov, nullptr);
ASSERT_NE(condAsimov, nullptr);

EXPECT_EQ(uncondObs->status(), 0);
EXPECT_EQ(condObs->status(), 0);
EXPECT_EQ(uncondAsimov->status(), 0);
EXPECT_EQ(condAsimov->status(), 0);

// The best-fit POI value from the unconditional fit result must be
// consistent with the one stored by the calculator.
auto *sFitUncond = static_cast<RooRealVar *>(uncondObs->floatParsFinal().find("s"));
ASSERT_NE(sFitUncond, nullptr);
EXPECT_DOUBLE_EQ(sFitUncond->getVal(), calc.GetMuHat()->getVal());

// In the conditional fits the POI is fixed to the tested value from the
// null-model snapshot, so it appears in the constant parameter list.
auto *sFitCond = static_cast<RooRealVar *>(condObs->constPars().find("s"));
ASSERT_NE(sFitCond, nullptr);
EXPECT_DOUBLE_EQ(sFitCond->getVal(), 50.0);

// The profile likelihood ratio test statistic reconstructed from the stored
// fit results must be non-negative, up to the same numerical tolerance that
// the calculator itself uses for qmu.
EXPECT_GE(condObs->minNll() - uncondObs->minNll(), -1.e-3);
EXPECT_GE(condAsimov->minNll() - uncondAsimov->minNll(), -1.e-3);
}
Loading