Proper way to use createNLL

Dear RooStats experts,

Given a workspace containing a mixture of constrained and unconstrained nuisance parameters, what is the appropriate way to use the createNLL method? Mainly I am unsure about are whether to include the Constrain and GlobalObservables options.

In the community, there seems to be a few common choices:

ModelConfig *mc = ...;
RooAbsData *data = ...;
RooAbsPdf *pdf=mc->GetPdf();
# choice 1
RooAbsReal* nll=pdf->createNLL(*data, 
				 Constrain(*mc->GetNuisanceParameters()),
				 GlobalObservables(*mc->GetGlobalObservables()), ...);
# choice 2
RooAbsReal* nll=pdf->createNLL(*data, 
				 GlobalObservables(*mc->GetGlobalObservables()), ...);

# choice 3
RooAbsReal* nll=pdf->createNLL(*data, 
				 Constrain(*mc->GetNuisanceParameters()));

I have tested the above choices and there are very slight differences in the resulting NLL value. From my understanding, *mc->GetNuisanceParameters() will return the set of nuisance parameters in the workspace whether or not they have the corresponding constraint pdfs. However only the constrained nuisance parameters have a corresponding global observable. I have no idea how createNLL handles this internally so I want to make sure I am not doing things incorrectly.

Welcome to the ROOT Forum!
I think @moneta can help

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.

Hi @AlkaidCheng!

Sorry for the late reply. Before writing more about your specific question, here are some good references about constraint terms and global observables:

  1. The documentation of RooAbsPdf::createNLL(), especially the explanation for the Constran() and GlobalObservables() command arguments
  2. An excellent forum post by @stephanh about global observables and constrained parameters
  3. For experts: the source code of the createConstraintTerm function that is also used by createNLL to double check what really happens
  4. RooFit tutorial on parameter constraints

Now, what is the right choice for you? Let’s look at your propositions again, noting what are used as default constraint parameters and global observables if you don’t pass them explicitly:

ModelConfig *mc = ...;
RooAbsData *data = ...;
RooAbsPdf *pdf=mc->GetPdf();
std::unique_ptr<RooArgSet> allParams{pdf->getParameters(*data)};

# choice 1
RooAbsReal* nll=pdf->createNLL(*data, 
				 Constrain(*mc->GetNuisanceParameters()),
				 GlobalObservables(*mc->GetGlobalObservables()),
               ...);
# choice 2
RooAbsReal* nll=pdf->createNLL(*data,
               // Constrained parameters not specified, so all parameters
               // are constrained if there are constraint terms for them.
               // It's like:
               // Constrain(*allParams),
				 GlobalObservables(*mc->GetGlobalObservables()),
               ...);

# choice 3
RooAbsReal* nll=pdf->createNLL(*data, 
				 Constrain(*mc->GetNuisanceParameters()),
               // GlobalObservables not specified,
               // so the constrained parameters are used, i.e. it's like:
               // GlobalObservables(*mc->GetNuisanceParameters()),
               ...);

Remember that the global observables are the variables that the constraint term is integrated over to determine its normalization. So not setting these global observables correctly is wrong if your constraints are not symmetric under the exchange of global observables and nuisance parameters! However, most constraints in the wild are Gaussian, so often you don’t notice that mistake, which is why you sill see your choice 3 in the wild even though it is conceptually wrong.

Choice 1 and 2 are both correct. In one you are just passing the set of constrained parameters explicitly, in the other not. Passing them explicitly can be useful if you study how removing some of the constraints will affect your fit result.

Just for you to try out, I have created also an example that showcases how you get a different fit result in a model with a Poisson constraint if you forget to set the global observables correctly:

void constraintsOptionsDemo()
{
   using namespace RooFit;

   // Construct a Gaussian pdf
   RooRealVar x("x", "x", -10, 10);
   RooRealVar m("m", "m", 5, 0, 10);
   RooRealVar s("s", "s", 2, 0.1, 10);
   RooGaussian gauss("gauss", "gauss(x,m,s)", x, m, s);

   // To reset parameters after fit to not bias results of the next fit
   auto resetParameters = [&](){
       m.setVal(5);
       m.setError(0);
       s.setVal(2);
       s.setError(0);
   };
 
   // Generate small dataset for our tests
   std::unique_ptr<RooDataSet> d{gauss.generate(x, 50)};

   // Create the global observables, which is constant (not a parameter)
   RooRealVar mGlobalObs("mGlobalObs", "mGlobalObs", 4, 0, 10);
   mGlobalObs.setConstant();
 
   // Construct Gaussian constraint pdf on parameter m
   RooPoisson constraint("constraint", "constraint", mGlobalObs, m);

   // Alternatively, you can try a Gaussian constraint, and you will see that
   // not defining the correct global observable doesn't matter here because
   // the integral of a Gaussian is symmetric under exchange of x and my.
   // RooGaussian constraint("constraint", "constraint", mGlobalObs, m, RooConst(2));
 
   // Multiply constraint with pdf
   RooProdPdf model("modec", "model with constraint", RooArgSet(gauss, constraint));
 
   // Create NLLs with different options
   std::unique_ptr<RooAbsReal> nll1{model.createNLL(*d)};
   std::unique_ptr<RooAbsReal> nll2{model.createNLL(*d, Constrain(m))};
   std::unique_ptr<RooAbsReal> nll3{model.createNLL(*d, Constrain(m), GlobalObservables(mGlobalObs))};

   // Do some fitting now, to show which differences in the NLL are not just an
   // absolute offset with no influence on the fit result.
   using R = std::unique_ptr<RooFitResult>;
   R r1{model.fitTo(*d, Save(), PrintLevel(-1))};
   resetParameters();
   R r2{model.fitTo(*d, Constrain(m), Save(), PrintLevel(-1))};
   resetParameters();
   R r3{model.fitTo(*d, Constrain(m), GlobalObservables(mGlobalObs), Save(), PrintLevel(-1))};
   resetParameters();

   // Print the likelihood values
   std::cout << nll1->getVal() << std::endl;
   std::cout << nll2->getVal() << std::endl;
   std::cout << nll3->getVal() << std::endl;

 
   // Print the fit results
   r1->Print();
   r2->Print();
   r3->Print();
}

Here is the output:

103.932
103.932
103.948

  RooFitResult: minimized FCN value: 103.81, estimated distance to minimum: 9.02216e-08
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 HESSE=0 

    Floating Parameter    FinalValue +/-  Error   
  --------------------  --------------------------
                     m    5.0597e+00 +/-  2.74e-01
                     s    1.9075e+00 +/-  2.06e-01


  RooFitResult: minimized FCN value: 103.81, estimated distance to minimum: 8.64109e-08
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 HESSE=0 

    Floating Parameter    FinalValue +/-  Error   
  --------------------  --------------------------
                     m    5.0597e+00 +/-  2.74e-01
                     s    1.9075e+00 +/-  2.06e-01


  RooFitResult: minimized FCN value: 103.825, estimated distance to minimum: 8.50857e-08
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 HESSE=0 

    Floating Parameter    FinalValue +/-  Error   
  --------------------  --------------------------
                     m    5.0612e+00 +/-  2.74e-01
                     s    1.9076e+00 +/-  2.06e-01

I hope this clarifies a few things!

Cheers,
Jonas