Dear Rooters,
I wanted to use the AsymptoticError flag in RooAbsPdf::fitTo without parameter limits, as Minuit transforms the parameters. The Minuit Manual states:
Because this transformation is non-linear, it is recommended
to avoid putting limits on parameters where they are not needed.
Parameter limits can lead to wrongly estimated errors and increased additional numerical inaccuracies (s. sec 1.3.1 and 5.3.2 in the manual linked above)
The problem can reproduced with the example provided by Christoph Langenbruch simply by dropping the range limits.
Here is the slightly modified version: rf611_weightedfits_AsymErrRng.C (9.8 KB)
Run with .x rf611_weightedfits_AsymErrRng.C(2, kFALSE) in the Root prompt. The second parameter is disabling the range.
When running without parameter limits one gets tons of warnings like:
[#0] ERROR:Eval -- RooAbsReal::logEvalError(pol) evaluation error,
origin : RooPolynomial::pol[ x=costheta coefList=(c0,c1) ]
message : p.d.f value is less than zero (-773502879776060572245688320.000000), forcing value to zero
server values: x=costheta=0.386751, coefList=(c0 = -2e+27 +/- 0.0183967,c1 = 0.0528082 +/- 0.036377)
Note the crazy value of -2e+27 for c0.
As I see it, this does not happen in the fit, but somehow in the calculation of the corrected error matrix.
So setting a small step size for Minuit would not help, as the parameters are not running away in the fit.
I am fearing that the parameter limits could lead to wrongly estimated errors by Minuit. Therefore also the asymptotically correct errors would be wrong, as they got a wrong input.
Is there a good reason why limits have to set for parameters when using AsymptoticError flag?
Setup
- Root 6.20.04 compiled with gcc 7.5
- Platform: Ubuntu 18.04
Appendix
I can not find any reason in the code I consider relevant. For convenience I quote the lines 1695-1746 from RooAbsPdf.cxx of tag v6-20-04:
if (doAsymptotic==1 && m.getNPar()>0) {
//Calculated corrected errors for weighted likelihood fits
std::unique_ptr<RooFitResult> rw(m.save());
//Weighted inverse Hessian matrix
const TMatrixDSym& matV = rw->covarianceMatrix();
coutI(Fitting) << "RooAbsPdf::fitTo(" << GetName() << ") Calculating covariance matrix according to the asymptotically correct approach. If you find this method useful please consider citing https://arxiv.org/abs/1911.01303." << endl;
//Initialise matrix containing first derivatives
TMatrixDSym num(rw->floatParsFinal().getSize());
for (int k=0; k<rw->floatParsFinal().getSize(); k++)
for (int l=0; l<rw->floatParsFinal().getSize(); l++)
num(k,l) = 0.0;
RooArgSet* obs = getObservables(data);
//Create derivative objects
std::vector<std::unique_ptr<RooDerivative> > derivatives;
const RooArgList& floated = rw->floatParsFinal();
std::unique_ptr<RooArgSet> floatingparams( (RooArgSet*)getParameters(data)->selectByAttrib("Constant", false) );
for (int k=0; k<floated.getSize(); k++) {
RooRealVar* paramresult = (RooRealVar*)floated.at(k);
RooRealVar* paraminternal = (RooRealVar*)floatingparams->find(paramresult->getTitle());
std::unique_ptr<RooDerivative> deriv( derivative(*paraminternal, *obs, 1) );
derivatives.push_back(std::move(deriv));
}
//Loop over data
for (int j=0; j<data.numEntries(); j++) {
//Sets obs to current data point, this is where the pdf will be evaluated
*obs = *data.get(j);
//Determine first derivatives
std::vector<Double_t> diffs(floated.getSize(), 0.0);
for (int k=0; k<floated.getSize(); k++) {
RooRealVar* paramresult = (RooRealVar*)floated.at(k);
RooRealVar* paraminternal = (RooRealVar*)floatingparams->find(paramresult->getTitle());
//First derivative to parameter k at best estimate point for this measurement
Double_t diff = derivatives.at(k)->getVal();
//Need to reset to best fit point after differentiation
*paraminternal = paramresult->getVal();
diffs.at(k) = diff;
}
//Fill numerator matrix
Double_t prob = getVal(obs);
for (int k=0; k<floated.getSize(); k++) {
for (int l=0; l<floated.getSize(); l++) {
num(k,l) += data.weight()*data.weight()*diffs.at(k)*diffs.at(l)/(prob*prob);
}
}
}
num.Similarity(matV);
//Propagate corrected errors to parameters objects
m.applyCovarianceMatrix(num);
}
