Likelihood and minimization: Recover from regions where the function is not defined.
We demonstrate improved recovery from disallowed parameters. For this, we use a polynomial PDF of the form \mathrm{Pol2} = \mathcal{N} \left( c + a_1 \cdot x + a_2 \cdot x^2 + 0.01 \cdot x^3 \right), where \f$ \mathcal{N} \f$ is a normalisation factor. Unless the parameters are chosen carefully, this function can be negative, and hence, it cannot be used as a PDF. In this case, RooFit passes an error to the minimiser, which might try to recover.
Author: Harshal Shende, Stephan Hageboeck (C++ version)
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Tuesday, March 19, 2024 at 07:17 PM.
import ROOT
Create a fit model: The polynomial is notoriously unstable, because it can quickly go negative. Since PDFs need to be positive, one often ends up with an unstable fit model.
x = ROOT.RooRealVar("x", "x", -15, 15)
a1 = ROOT.RooRealVar("a1", "a1", -0.5, -10.0, 20.0)
a2 = ROOT.RooRealVar("a2", "a2", 0.2, -10.0, 20.0)
a3 = ROOT.RooRealVar("a3", "a3", 0.01)
pdf = ROOT.RooPolynomial("pol3", "c + a1 * x + a2 * x*x + 0.01 * x*x*x", x, [a1, a2, a3])
Create toy data with all-positive coefficients:
data = pdf.generate(x, 10000)
For plotting. We create pointers to the plotted objects. We want these objects to leak out of the function, so we can still see them after it returns.
c = ROOT.TCanvas()
frame = x.frame()
data.plotOn(frame, Name="data")
<cppyy.gbl.RooPlot object at 0xaee0440>
Plotting a PDF with disallowed parameters doesn't work. We would get a lot of error messages. Therefore, we disable plotting messages in RooFit's message streams:
ROOT.RooMsgService.instance().getStream(0).removeTopic(ROOT.RooFit.Plotting)
ROOT.RooMsgService.instance().getStream(1).removeTopic(ROOT.RooFit.Plotting)
Before 6.24, RooFit wasn't able to recover from invalid parameters. The minimiser just errs around the starting values of the parameters without finding any improvement.
Set up the parameters such that the PDF would come out negative. The PDF is now undefined.
a1.setVal(10.0)
a2.setVal(-1.0)
Perform a fit:
fitWithoutRecovery = pdf.fitTo(
data,
Save=True,
RecoverFromUndefinedRegions=0.0, # This is how RooFit behaved prior to ROOT 6.24
PrintEvalErrors=-1, # We are expecting a lot of evaluation errors. -1 switches off printing.
PrintLevel=-1,
)
pdf.plotOn(frame, LineColor="r", Name="noRecovery")
<cppyy.gbl.RooPlot object at 0xaee0440>
[#1] INFO:Fitting -- RooAbsPdf::fitTo(pol3_over_pol3_Int[x]) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- using CPU computation library compiled with -mavx2 [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_pol3_over_pol3_Int[x]_pol3Data) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization [#0] ERROR:Eval -- RooAbsReal::logEvalError(pol3) evaluation error, origin : RooPolynomial::pol3[ x=x coefList=(a1,a2,a3) ] message : p.d.f normalization integral is zero or negative: -2220.000000 server values: x=x=0, coefList=(a1 = 10 +/- 0,a2 = -1 +/- 0,a3 = 0.01)
Warning in <ROOT::Math::Fitter::CalculateHessErrors>: Error when calculating Hessian
The minimiser gets information about the "badness" of the violation of the function definition. It uses this to find its way out of the disallowed parameter regions.
print("\n\n\n-------------- Starting second fit ---------------\n\n")
-------------- Starting second fit ---------------
Reset the parameters such that the PDF is again undefined.
a1.setVal(10.0)
a2.setVal(-1.0)
Fit again, but pass recovery information to the minimiser:
fitWithRecovery = pdf.fitTo(
data,
Save=True,
RecoverFromUndefinedRegions=1.0, # The magnitude of the recovery information can be chosen here.
# Higher values mean more aggressive recovery.
PrintEvalErrors=-1, # We are still expecting a few evaluation errors.
PrintLevel=0,
)
pdf.plotOn(frame, LineColor="b", Name="recovery")
<cppyy.gbl.RooPlot object at 0xaee0440>
[#1] INFO:Fitting -- RooAbsPdf::fitTo(pol3_over_pol3_Int[x]) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_pol3_over_pol3_Int[x]_pol3Data) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization Minuit2Minimizer: Minimize with max-calls 1000 convergence for edm < 1 strategy 1 Minuit2Minimizer : Valid minimum - status = 0 FVAL = -1002.2262595660759 Edm = 2.95538313214564806e-09 Nfcn = 251 a1 = -0.498159 +/- 0.0227242 (limited) a2 = 0.198316 +/- 0.00564906 (limited) [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
Warning in <Minuit2>: MnPosDef non-positive diagonal element in covariance matrix[ 0 ] = -13.133 Warning in <Minuit2>: MnPosDef Added to diagonal of Error matrix a value 13.633 Warning in <Minuit2>: DavidonErrorUpdator delgam < 0 : first derivatives increasing along search line (details in info log) Warning in <Minuit2>: VariableMetricBuilder Matrix not pos.def, gdel = 67.8685 > 0 Warning in <Minuit2>: MnPosDef non-positive diagonal element in covariance matrix[ 0 ] = -0.000525791 Warning in <Minuit2>: MnPosDef non-positive diagonal element in covariance matrix[ 1 ] = -0.00482974 Warning in <Minuit2>: MnPosDef Added to diagonal of Error matrix a value 0.50483 Warning in <Minuit2>: VariableMetricBuilder gdel = -51905.4 Warning in <Minuit2>: DavidonErrorUpdator delgam < 0 : first derivatives increasing along search line (details in info log) Warning in <Minuit2>: VariableMetricBuilder Matrix not pos.def, gdel = 56.5625 > 0 Warning in <Minuit2>: MnPosDef non-positive diagonal element in covariance matrix[ 1 ] = -0.00239309 Warning in <Minuit2>: MnPosDef Added to diagonal of Error matrix a value 0.502393 Warning in <Minuit2>: VariableMetricBuilder gdel = -15405 Warning in <Minuit2>: DavidonErrorUpdator delgam < 0 : first derivatives increasing along search line (details in info log) Warning in <Minuit2>: VariableMetricBuilder Matrix not pos.def, gdel = 3573.84 > 0 Warning in <Minuit2>: MnPosDef non-positive diagonal element in covariance matrix[ 0 ] = -1.28847e-05 Warning in <Minuit2>: MnPosDef non-positive diagonal element in covariance matrix[ 1 ] = -2.85488e-05 Warning in <Minuit2>: MnPosDef Added to diagonal of Error matrix a value 0.500029 Warning in <Minuit2>: VariableMetricBuilder gdel = -5.37564e+07
We print the two fit results, and plot the fitted curves. The curve of the fit without recovery cannot be plotted, because the PDF is undefined if a2 < 0.
fitWithoutRecovery.Print()
print(
"Without recovery, the fitter encountered {}".format(fitWithoutRecovery.numInvalidNLL())
+ " invalid function values. The parameters are unchanged.\n"
)
fitWithRecovery.Print()
print(
"With recovery, the fitter encountered {}".format(fitWithoutRecovery.numInvalidNLL())
+ " invalid function values, but the parameters are fitted.\n"
)
legend = ROOT.TLegend(0.5, 0.7, 0.9, 0.9)
legend.SetBorderSize(0)
legend.SetFillStyle(0)
legend.AddEntry("data", "Data", "P")
legend.AddEntry("noRecovery", "Without recovery (cannot be plotted)", "L")
legend.AddEntry("recovery", "With recovery", "L")
frame.Draw()
legend.Draw()
c.Draw()
c.SaveAs("rf612_recoverFromInvalidParameters.png")
Without recovery, the fitter encountered 23 invalid function values. The parameters are unchanged. With recovery, the fitter encountered 23 invalid function values, but the parameters are fitted. RooFitResult: minimized FCN value: 0, estimated distance to minimum: 0 covariance matrix quality: Not calculated at all Status : MINIMIZE=-1 HESSE=302 Floating Parameter FinalValue +/- Error -------------------- -------------------------- a1 1.0000e+01 +/- 0.00e+00 a2 -1.0000e+00 +/- 0.00e+00 RooFitResult: minimized FCN value: 29650.9, estimated distance to minimum: 2.95925e-09 covariance matrix quality: Full, accurate covariance matrix Status : MINIMIZE=0 HESSE=0 Floating Parameter FinalValue +/- Error -------------------- -------------------------- a1 -4.9816e-01 +/- 2.27e-02 a2 1.9832e-01 +/- 5.65e-03
Info in <TCanvas::Print>: png file rf612_recoverFromInvalidParameters.png has been created
Draw all canvases
from ROOT import gROOT
gROOT.GetListOfCanvases().Draw()