|
|
|
#include <FitterAlgoBase.h>
Public Member Functions | |
| void | applyOptionsBase (const boost::program_options::variables_map &vm) |
| FitterAlgoBase (const char *title="<FillMe> specific options") | |
| virtual bool | run (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint) |
Protected Member Functions | |
| RooFitResult * | doFit (RooAbsPdf &pdf, RooAbsData &data, RooRealVar &r, const RooCmdArg &constrain, bool doHesse=true, int ndim=1, bool reuseNLL=false) |
| RooFitResult * | doFit (RooAbsPdf &pdf, RooAbsData &data, const RooArgList &rs, const RooCmdArg &constrain, bool doHesse=true, int ndim=1, bool reuseNLL=false) |
| double | findCrossing (CascadeMinimizer &minim, RooAbsReal &nll, RooRealVar &r, double level, double rStart, double rBound) |
| virtual bool | runSpecific (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint)=0 |
Protected Attributes | |
| static bool | do95_ = false |
| static bool | keepFailures_ = false |
| static std::string | minimizerAlgoForMinos_ = "Minuit2,simplex" |
| static int | minimizerStrategyForMinos_ = 0 |
| static float | minimizerToleranceForMinos_ = 1e-4 |
| std::auto_ptr< RooAbsReal > | nll |
Static Protected Attributes | |
| static int | maxFailedSteps_ = 5 |
| static std::string | minimizerAlgo_ = "Minuit2" |
| static int | minimizerStrategy_ = 1 |
| static float | minimizerTolerance_ = 1e-2 |
| static float | nllValue_ = std::numeric_limits<float>::quiet_NaN() |
| static float | preFitValue_ = 1.0 |
| static bool | robustFit_ = false |
| static bool | saveNLL_ = false |
| static float | stepSize_ = 0.1 |
Do a ML fit of the data with background and signal+background hypothesis and print out diagnostics plots
Definition at line 20 of file FitterAlgoBase.h.
| FitterAlgoBase::FitterAlgoBase | ( | const char * | title = "<FillMe> specific options" | ) |
Definition at line 53 of file FitterAlgoBase.cc.
References do95_, maxFailedSteps_, minimizerAlgo_, minimizerAlgoForMinos_, minimizerStrategy_, minimizerStrategyForMinos_, minimizerTolerance_, minimizerToleranceForMinos_, LimitAlgo::options_, preFitValue_, robustFit_, and stepSize_.
:
LimitAlgo(title)
{
options_.add_options()
("minimizerAlgo", boost::program_options::value<std::string>(&minimizerAlgo_)->default_value(minimizerAlgo_), "Choice of minimizer (Minuit vs Minuit2)")
("minimizerTolerance", boost::program_options::value<float>(&minimizerTolerance_)->default_value(minimizerTolerance_), "Tolerance for minimizer")
("minimizerStrategy", boost::program_options::value<int>(&minimizerStrategy_)->default_value(minimizerStrategy_), "Stragegy for minimizer")
("preFitValue", boost::program_options::value<float>(&preFitValue_)->default_value(preFitValue_), "Value of signal strength pre-fit")
("do95", boost::program_options::value<bool>(&do95_)->default_value(do95_), "Compute also 2-sigma interval from delta(nll) = 1.92 instead of 0.5")
("robustFit", boost::program_options::value<bool>(&robustFit_)->default_value(robustFit_), "Search manually for 1 and 2 sigma bands instead of using Minos")
("maxFailedSteps", boost::program_options::value<int>(&maxFailedSteps_)->default_value(maxFailedSteps_), "How many failed steps to retry before giving up")
("stepSize", boost::program_options::value<float>(&stepSize_)->default_value(stepSize_), "Step size for robust fits (multiplier of the range)")
("minimizerAlgoForMinos", boost::program_options::value<std::string>(&minimizerAlgoForMinos_)->default_value(minimizerAlgoForMinos_), "Choice of minimizer (Minuit vs Minuit2) for profiling in robust fits")
("minimizerStrategyForMinos", boost::program_options::value<int>(&minimizerStrategyForMinos_)->default_value(minimizerStrategyForMinos_), "Stragegy for minimizer for profiling in robust fits")
("minimizerToleranceForMinos", boost::program_options::value<float>(&minimizerToleranceForMinos_)->default_value(minimizerToleranceForMinos_), "Tolerance for minimizer for profiling in robust fits")
("saveNLL", "Save the negative log-likelihood at the minimum in the output tree (note: value is relative to the pre-fit state)")
("keepFailures", "Save the results even if the fit is declared as failed (for NLL studies)")
;
}
| void FitterAlgoBase::applyOptionsBase | ( | const boost::program_options::variables_map & | vm | ) |
Definition at line 73 of file FitterAlgoBase.cc.
References keepFailures_, and saveNLL_.
Referenced by ChannelCompatibilityCheck::applyOptions(), MultiDimFit::applyOptions(), and MaxLikelihoodFit::applyOptions().
{
saveNLL_ = vm.count("saveNLL");
keepFailures_ = vm.count("keepFailures");
}
| RooFitResult * FitterAlgoBase::doFit | ( | RooAbsPdf & | pdf, |
| RooAbsData & | data, | ||
| RooRealVar & | r, | ||
| const RooCmdArg & | constrain, | ||
| bool | doHesse = true, |
||
| int | ndim = 1, |
||
| bool | reuseNLL = false |
||
| ) | [protected] |
Fit data with pdf, with parameters of interest in r, and specified constraint If ndim = 1, errors on each parameter are from a 1-dim chisquare, as for a single parameter fit If ndim > 1, errors on each parameter are from a n-dim chisquare, as for a joint estimation of N parameters
Definition at line 90 of file FitterAlgoBase.cc.
Referenced by ChannelCompatibilityCheck::runSpecific(), MultiDimFit::runSpecific(), and MaxLikelihoodFit::runSpecific().
| RooFitResult * FitterAlgoBase::doFit | ( | RooAbsPdf & | pdf, |
| RooAbsData & | data, | ||
| const RooArgList & | rs, | ||
| const RooCmdArg & | constrain, | ||
| bool | doHesse = true, |
||
| int | ndim = 1, |
||
| bool | reuseNLL = false |
||
| ) | [protected] |
Definition at line 94 of file FitterAlgoBase.cc.
References CloseCoutSentry::clear(), CascadeMinimizer::Constrained, gather_cfg::cout, do95_, findCrossing(), i, CascadeMinimizer::improve(), edm::detail::isnan(), keepFailures_, max(), CascadeMinimizer::minimize(), CascadeMinimizer::minimizer(), minimizerStrategy_, minimizerStrategyForMinos_, n, nll, nllValue_, convertSQLiteXML::ok, alignCSCRings::r, run_regression::ret, robustFit_, CascadeMinimizer::save(), CascadeMinimizer::setErrorLevel(), CascadeMinimizer::setStrategy(), and CascadeMinimizer::Unconstrained.
{
RooFitResult *ret = 0;
if (reuseNLL) nll->setData(data); // reuse nll but swap out the data
else nll.reset(pdf.createNLL(data, constrain, RooFit::Extended(pdf.canBeExtended()))); // make a new nll
double nll0 = nll->getVal();
double delta68 = 0.5*ROOT::Math::chisquared_quantile_c(1-0.68,ndim);
double delta95 = 0.5*ROOT::Math::chisquared_quantile_c(1-0.95,ndim);
CascadeMinimizer minim(*nll, CascadeMinimizer::Unconstrained, rs.getSize() ? dynamic_cast<RooRealVar*>(rs.first()) : 0);
minim.setStrategy(minimizerStrategy_);
minim.setErrorLevel(delta68);
CloseCoutSentry sentry(verbose < 3);
bool ok = minim.minimize(verbose-1);
nllValue_ = nll->getVal() - nll0;
if (!ok && !keepFailures_) { std::cout << "Initial minimization failed. Aborting." << std::endl; return 0; }
if (doHesse) minim.minimizer().hesse();
sentry.clear();
ret = minim.save();
if (verbose > 1) { ret->Print("V"); }
std::auto_ptr<RooArgSet> allpars(pdf.getParameters(data));
for (int i = 0, n = rs.getSize(); i < n; ++i) {
// if this is not the first fit, reset parameters
if (i) {
RooArgSet oldparams(ret->floatParsFinal());
*allpars = oldparams;
}
// get the parameter to scan, amd output variable in fit result
RooRealVar &r = dynamic_cast<RooRealVar &>(*rs.at(i));
RooRealVar &rf = dynamic_cast<RooRealVar &>(*ret->floatParsFinal().find(r.GetName()));
double r0 = r.getVal(), rMin = r.getMin(), rMax = r.getMax();
if (!robustFit_) {
if (do95_) {
throw std::runtime_error("95% CL errors with Minos are not working at the moment.");
minim.setErrorLevel(delta95);
minim.improve(verbose-1);
minim.setErrorLevel(delta95);
if (minim.minimizer().minos(RooArgSet(r)) != -1) {
rf.setRange("err95", r.getVal() + r.getAsymErrorLo(), r.getVal() + r.getAsymErrorHi());
}
minim.setErrorLevel(delta68);
minim.improve(verbose-1);
}
if (minim.minimizer().minos(RooArgSet(r)) != -1) {
rf.setRange("err68", r.getVal() + r.getAsymErrorLo(), r.getVal() + r.getAsymErrorHi());
rf.setAsymError(r.getAsymErrorLo(), r.getAsymErrorHi());
}
} else {
r.setVal(r0); r.setConstant(true);
CascadeMinimizer minim2(*nll, CascadeMinimizer::Constrained);
minim2.setStrategy(minimizerStrategyForMinos_);
std::auto_ptr<RooArgSet> allpars(nll->getParameters((const RooArgSet *)0));
double nll0 = nll->getVal();
double threshold68 = nll0 + delta68;
double threshold95 = nll0 + delta95;
// search for crossings
assert(!std::isnan(r0));
// high error
double hi68 = findCrossing(minim2, *nll, r, threshold68, r0, rMax);
double hi95 = do95_ ? findCrossing(minim2, *nll, r, threshold95, std::isnan(hi68) ? r0 : hi68, std::max(rMax, std::isnan(hi68*2-r0) ? r0 : hi68*2-r0)) : r0;
// low error
*allpars = RooArgSet(ret->floatParsFinal()); r.setVal(r0); r.setConstant(true);
double lo68 = findCrossing(minim2, *nll, r, threshold68, r0, rMin);
double lo95 = do95_ ? findCrossing(minim2, *nll, r, threshold95, std::isnan(lo68) ? r0 : lo68, rMin) : r0;
rf.setAsymError(!std::isnan(lo68) ? lo68 - r0 : 0, !std::isnan(hi68) ? hi68 - r0 : 0);
rf.setRange("err68", !std::isnan(lo68) ? lo68 : r0, !std::isnan(hi68) ? hi68 : r0);
if (do95_ && (!std::isnan(lo95) || !std::isnan(hi95))) {
rf.setRange("err95", !std::isnan(lo95) ? lo95 : r0, !std::isnan(hi95) ? hi95 : r0);
}
r.setVal(r0); r.setConstant(false);
}
}
return ret;
}
| double FitterAlgoBase::findCrossing | ( | CascadeMinimizer & | minim, |
| RooAbsReal & | nll, | ||
| RooRealVar & | r, | ||
| double | level, | ||
| double | rStart, | ||
| double | rBound | ||
| ) | [protected] |
Definition at line 179 of file FitterAlgoBase.cc.
References gather_cfg::cout, runtimedef::get(), CascadeMinimizer::improve(), max(), maxFailedSteps_, min, minimizerAlgoForMinos_, minimizerToleranceForMinos_, convertSQLiteXML::ok, CascadeMinimizer::save(), and stepSize_.
Referenced by MultiDimFit::doBox(), MultiDimFit::doContour2D(), and doFit().
{
ProfileLikelihood::MinimizerSentry minimizerConfig(minimizerAlgoForMinos_, minimizerToleranceForMinos_);
if (verbose) std::cout << "Searching for crossing at nll = " << level << " in the interval " << rStart << ", " << rBound << std::endl;
double rInc = stepSize_*(rBound - rStart);
r.setVal(rStart);
std::auto_ptr<RooFitResult> checkpoint;
std::auto_ptr<RooArgSet> allpars;
bool ok = false;
{
CloseCoutSentry sentry(verbose < 3);
ok = minim.improve(verbose-1);
checkpoint.reset(minim.save());
}
if (!ok) { std::cout << "Error: minimization failed at " << r.GetName() << " = " << rStart << std::endl; return NAN; }
double here = nll.getVal();
int nfail = 0;
if (verbose > 0) { printf(" %s lvl-here lvl-there stepping\n", r.GetName()); fflush(stdout); }
do {
rStart += rInc;
if (rInc*(rStart - rBound) > 0) { // went beyond bounds
rStart -= rInc;
rInc = 0.5*(rBound-rStart);
}
r.setVal(rStart);
nll.clearEvalErrorLog(); nll.getVal();
if (nll.numEvalErrors() > 0) {
ok = false;
} else {
CloseCoutSentry sentry(verbose < 3);
ok = minim.improve(verbose-1);
}
if (!ok) {
nfail++;
if (nfail >= maxFailedSteps_) { std::cout << "Error: minimization failed at " << r.GetName() << " = " << rStart << std::endl; return NAN; }
RooArgSet oldparams(checkpoint->floatParsFinal());
if (allpars.get() == 0) allpars.reset(nll.getParameters((const RooArgSet *)0));
*allpars = oldparams;
rStart -= rInc; rInc *= 0.5;
continue;
} else nfail = 0;
double there = here;
here = nll.getVal();
if (verbose > 0) { printf("%f %+.5f %+.5f %f\n", rStart, level-here, level-there, rInc); fflush(stdout); }
if ( fabs(here - level) < 4*minimizerToleranceForMinos_ ) {
// set to the right point with interpolation
r.setVal(rStart + (level-here)*(level-there)/(here-there));
return r.getVal();
} else if (here > level) {
// I'm above the level that I wanted, this means I stepped too long
// First I take back all the step
rStart -= rInc;
// Then I try to figure out a better step
if (runtimedef::get("FITTER_DYN_STEP")) {
if (fabs(there - level) > 0.05) { // If started from far away, I still have to step carefully
double rIncFactor = std::max(0.2, std::min(0.7, 0.75*(level-there)/(here-there)));
//printf("\t\t\t\t\tCase A1: level-there = %f, here-there = %f, rInc(Old) = %f, rInFactor = %f, rInc(New) = %f\n", level-there, here-there, rInc, rIncFactor, rInc*rIncFactor);
rInc *= rIncFactor;
} else { // close enough to jump straight to target
double rIncFactor = std::max(0.05, std::min(0.95, 0.95*(level-there)/(here-there)));
//printf("\t\t\t\t\tCase A2: level-there = %f, here-there = %f, rInc(Old) = %f, rInFactor = %f, rInc(New) = %f\n", level-there, here-there, rInc, rIncFactor, rInc*rIncFactor);
rInc *= rIncFactor;
}
} else {
rInc *= 0.3;
}
if (allpars.get() == 0) allpars.reset(nll.getParameters((const RooArgSet *)0));
RooArgSet oldparams(checkpoint->floatParsFinal());
*allpars = oldparams;
} else if ((here-there)*(level-there) < 0 && // went wrong
fabs(here-there) > 0.1) { // by more than roundoff
if (allpars.get() == 0) allpars.reset(nll.getParameters((const RooArgSet *)0));
RooArgSet oldparams(checkpoint->floatParsFinal());
*allpars = oldparams;
rStart -= rInc; rInc *= 0.5;
} else {
// I did this step, and I'm not there yet
if (runtimedef::get("FITTER_DYN_STEP")) {
if (fabs(here - level) > 0.05) { // we still have to step carefully
if ((here-there)*(level-there) > 0) { // if we went in the right direction
// then optimize step size
double rIncFactor = std::max(0.2, std::min(2.0, 0.75*(level-there)/(here-there)));
//printf("\t\t\t\t\tCase B1: level-there = %f, here-there = %f, rInc(Old) = %f, rInFactor = %f, rInc(New) = %f\n", level-there, here-there, rInc, rIncFactor, rInc*rIncFactor);
rInc *= rIncFactor;
} //else printf("\t\t\t\t\tCase B3: level-there = %f, here-there = %f, rInc(Old) = %f\n", level-there, here-there, rInc);
} else { // close enough to jump straight to target
double rIncFactor = std::max(0.05, std::min(4.0, 0.95*(level-there)/(here-there)));
//printf("\t\t\t\t\tCase B2: level-there = %f, here-there = %f, rInc(Old) = %f, rInFactor = %f, rInc(New) = %f\n", level-there, here-there, rInc, rIncFactor, rInc*rIncFactor);
rInc *= rIncFactor;
}
} else {
//nothing?
}
checkpoint.reset(minim.save());
}
} while (fabs(rInc) > minimizerToleranceForMinos_*stepSize_*std::max(1.0,rBound-rStart));
if (fabs(here - level) > 0.01) {
std::cout << "Error: closed range without finding crossing." << std::endl;
return NAN;
} else {
return r.getVal();
}
}
| bool FitterAlgoBase::run | ( | RooWorkspace * | w, |
| RooStats::ModelConfig * | mc_s, | ||
| RooStats::ModelConfig * | mc_b, | ||
| RooAbsData & | data, | ||
| double & | limit, | ||
| double & | limitErr, | ||
| const double * | hint | ||
| ) | [virtual] |
Implements LimitAlgo.
Definition at line 79 of file FitterAlgoBase.cc.
References Combine::addBranch(), minimizerAlgo_, minimizerTolerance_, nllValue_, runSpecific(), and saveNLL_.
{
ProfileLikelihood::MinimizerSentry minimizerConfig(minimizerAlgo_, minimizerTolerance_);
CloseCoutSentry sentry(verbose < 0);
static bool shouldCreateNLLBranch = saveNLL_;
if (shouldCreateNLLBranch) { Combine::addBranch("nll", &nllValue_, "nll/F"); shouldCreateNLLBranch = false; }
return runSpecific(w, mc_s, mc_b, data, limit, limitErr, hint);
}
| virtual bool FitterAlgoBase::runSpecific | ( | RooWorkspace * | w, |
| RooStats::ModelConfig * | mc_s, | ||
| RooStats::ModelConfig * | mc_b, | ||
| RooAbsData & | data, | ||
| double & | limit, | ||
| double & | limitErr, | ||
| const double * | hint | ||
| ) | [protected, pure virtual] |
Implemented in ChannelCompatibilityCheck, MaxLikelihoodFit, and MultiDimFit.
Referenced by run().
bool FitterAlgoBase::do95_ = false [protected] |
Definition at line 36 of file FitterAlgoBase.h.
Referenced by doFit(), MultiDimFit::doSingles(), FitterAlgoBase(), ChannelCompatibilityCheck::runSpecific(), and MaxLikelihoodFit::runSpecific().
bool FitterAlgoBase::keepFailures_ = false [protected] |
Definition at line 40 of file FitterAlgoBase.h.
Referenced by applyOptionsBase(), doFit(), and MultiDimFit::runSpecific().
int FitterAlgoBase::maxFailedSteps_ = 5 [static, protected] |
Definition at line 38 of file FitterAlgoBase.h.
Referenced by findCrossing(), and FitterAlgoBase().
std::string FitterAlgoBase::minimizerAlgo_ = "Minuit2" [static, protected] |
Definition at line 30 of file FitterAlgoBase.h.
Referenced by FitterAlgoBase(), and run().
std::string FitterAlgoBase::minimizerAlgoForMinos_ = "Minuit2,simplex" [protected] |
Definition at line 30 of file FitterAlgoBase.h.
Referenced by findCrossing(), and FitterAlgoBase().
int FitterAlgoBase::minimizerStrategy_ = 1 [static, protected] |
Definition at line 32 of file FitterAlgoBase.h.
Referenced by MultiDimFit::doBox(), MultiDimFit::doContour2D(), doFit(), MultiDimFit::doGrid(), MultiDimFit::doRandomPoints(), FitterAlgoBase(), and MaxLikelihoodFit::runSpecific().
int FitterAlgoBase::minimizerStrategyForMinos_ = 0 [protected] |
Definition at line 32 of file FitterAlgoBase.h.
Referenced by doFit(), and FitterAlgoBase().
float FitterAlgoBase::minimizerTolerance_ = 1e-2 [static, protected] |
Definition at line 31 of file FitterAlgoBase.h.
Referenced by FitterAlgoBase(), and run().
float FitterAlgoBase::minimizerToleranceForMinos_ = 1e-4 [protected] |
Definition at line 31 of file FitterAlgoBase.h.
Referenced by findCrossing(), and FitterAlgoBase().
std::auto_ptr<RooAbsReal> FitterAlgoBase::nll [protected] |
Definition at line 42 of file FitterAlgoBase.h.
Referenced by doFit(), MultiDimFit::runSpecific(), and MaxLikelihoodFit::runSpecific().
float FitterAlgoBase::nllValue_ = std::numeric_limits<float>::quiet_NaN() [static, protected] |
Definition at line 41 of file FitterAlgoBase.h.
float FitterAlgoBase::preFitValue_ = 1.0 [static, protected] |
Definition at line 34 of file FitterAlgoBase.h.
Referenced by FitterAlgoBase(), ChannelCompatibilityCheck::runSpecific(), and MaxLikelihoodFit::runSpecific().
bool FitterAlgoBase::robustFit_ = false [static, protected] |
Definition at line 36 of file FitterAlgoBase.h.
Referenced by doFit(), and FitterAlgoBase().
bool FitterAlgoBase::saveNLL_ = false [static, protected] |
Definition at line 40 of file FitterAlgoBase.h.
Referenced by applyOptionsBase(), and run().
float FitterAlgoBase::stepSize_ = 0.1 [static, protected] |
Definition at line 37 of file FitterAlgoBase.h.
Referenced by findCrossing(), and FitterAlgoBase().
1.7.3