RooStats hypothesis inversion demo, on/off problem, how to handle observables

Roofit and RooStats
1

Dear RooStats users,
I am quite new to the RooStats software, and I am trying to set a 90% CL upper limit on a signal strength parameter, that enters in an “on-off” problem.
My PDF reads:

P(n,n’;s,b) = Poisson(n,s+b) * Posson (n’,b)

where n and n’ are my observables, s is the POI, and b a nuisance parameter. I am simplifying the problem here, assuming that the ancillary measurement to determine b has the same luminosity / efficiency as the main one.

I have a doubt regarding how to handle n’. I see these two options:

OPTION-1: handle n’ as a “global fixed variable”

int nsb = 30;   // number of observed events
  int nb = 28; // number of background events
  
  RooWorkspace w("w",true);
  
  // make Poisson model
  w.factory("sum:nexp(s[2,0,20],b[0,0,50])");
  // Poisson of (nsb_obs | s+b)
  w.factory("Poisson:psb(nsb_obs[50,0,70],nexp)");
  // Poisson of (nb_obs | b)
  w.factory("Poisson:pb(nb_obs[50,0,50],b)");
  //The model
  w.factory("PROD:model(psb,pb)");
  
  
  w.var("nb_obs")->setVal(nb);
  w.var("nb_obs")->setConstant(true); // needed for being treated as global observables
 
  //set value of observed events
  RooRealVar *nsb_obs = w.var("nsb_obs");
  nsb_obs->setVal(nsb);

  // make data set with the number of observed events
  RooDataSet data("data","",*nsb_obs);
  data.add(*nsb_obs);
  w.import(data);

  w.Print();

  ModelConfig mc("Smodel",&w);
   mc.SetPdf(*w.pdf("model"));
   mc.SetParametersOfInterest(*w.var("s"));
   mc.SetObservables(*w.var("nsb_obs"));
  
   // these are needed for the hypothesis tests
   mc.SetSnapshot(*w.var("s"));
   mc.SetGlobalObservables(*w.var("nb_obs"));

OPTION-2: handle n’ as an observable

 // make Poisson model
  w.factory("sum:nexp(s[2,0,20],b[0,0,50])");
  // Poisson of (nsb_obs | s+b)
  w.factory("Poisson:psb(nsb_obs[50,0,70],nexp)");
  // Poisson of (nb_obs | b)
  w.factory("Poisson:pb(nb_obs[50,0,50],b)");
  //The model
  w.factory("PROD:model(psb,pb)");
   
  //set value of observed events
  RooRealVar *nsb_obs = w.var("nsb_obs");
  nsb_obs->setVal(nsb);

  //set value of observed events in ancillary measurement
  RooRealVar *nb_obs = w.var("nb_obs");
  nb_obs->setVal(nb);

  // make data set with the number of observed events in both measurements
  RooDataSet data("data","", RooArgSet(*nsb_obs,*nb_obs ));
  data.add(RooArgSet(*nsb_obs,*nb_obs ));
  w.import(data);

  w.Print();

   ModelConfig mc("Smodel",&w);
   mc.SetPdf(*w.pdf("model"));
   mc.SetParametersOfInterest(*w.var("s"));
   mc.SetObservables(RooArgSet(*w.var("nsb_obs"),*w.var("nb_obs")));

  mc.SetSnapshot(*w.var("s"));

With both these strategies, I can compute the 90% CL limit, using for example the StandardHypoTestInvDemo.C (I omitted the code to create the B-only model).

I get same results concerning the CL limit (CLs+b -> s=11.8464), but different values for the expected upper limits, using the B (alternate) model.

What is the difference between these two approaches? I think that, in this case, nb_obs should be declared as a “observables”, i.e. as a probabilistic variables?

Thanks,
Bests,

Andrea

2

I guess @StephanH can take a look and give you some details

3

I expect that it is throwing toys for n and n’ in the first case, and in the second case it’s only throwing toys for n.
I’m confused, though, why you register it as observable in the second case? It should also be a global observable. Global observables are things that you have to throw toys for, but they are not in a dataset. Normally, RooFit throws toys for everything that’s in a dataset, but not for other parameters. That’s why you call them “global” observable.

4

Hi StephanH,
thanks for your explanation.

However, why shouldn’t be n’ (I call it nb_obs, sorry for mis-naming) also an observable? I consider here two measurements, and I have both n and n’ in the dataset “data”, right? In other words, may I ask what is the difference between an observable and a “global observable”?

Thanks
Andrea

5

An observable is found in a dataset, so RooFit knows what to do with it.
A global observable is not in any dataset, so RooFit might mistake it for a parameter. That’s why we explicitly call it global observable. I don’t know how you throw the toys, but it could well be that it’s mistaking n’ for a parameter. Can you try to increase the print level or print the toy datasets?

I discussed with a colleague, and we agreed that “in principle” they should be equivalent.

6

I see now that I confused the two cases. You did exactly what I would have recommended, but I’m surprised that you get different results.

  • What are you using to throw the toys?
  • Could it be different random numbers? (ToyMCCalculator)

Maybe it’s good to post the full example.

7

Dear StephanH, sure. For simplicity, I post here the following codes.

  1. test2.C: the code that uses the case “1” (n’ as a global variable) to produce a TFile “test2.out1.root” with the workspace containing the S+B ModelConfig and the B ModelConfig

  2. test2a.C: the code that uses the case “2” (n’ as an observable) to produce a TFile “test2a.out1.root” with the workspace containing the S+B ModelConfig and the B ModelConfig

  3. test2HipoInverter.C: the code that performs the CL determination. At the moment, I am using a one-sided AsymptoticCalculator. There is also a commented part for a Frequentist Calculator. It is sufficient to change the name of the “infile” to switch between the two cases.

Thanks,
Bests
Andreatest2.C (1.6 KB) test2a.C (1.7 KB) test2HipoInverter.C (3.8 KB)

8

Hi,

Thank you for posting this issue. I have found a bug in handling the global observable in case of a Poisson term. The global observable is rounded wrongly to an integer in case of the Asimov fit. This explains the difference. The result with the two observable is instead correct.
I have opened a JIRA ticket:
https://sft.its.cern.ch/jira/browse/ROOT-10920

Best regards

Lorenzo

9

Hi,
This is now fixed in the ROOT master by this PR,


that is now merged in ROOT master

Lorenzo

10

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