As far as I understand, the ROOT macro StandardHypoTestInvDemo.C is specially designed to calculate a confidence interval for the particular case in which the NULL hypothesis is S+B (i.e. you are testing at which level of signal your model is incompatible with the data you observed).
I am trying to use the same code to test the B-only hypothesis, so that if if the null hypothesis is ruled out at a given confidence level, one could claim a discovery. In some tutorials online I’ve seen that in principle that should be as easy as swaping the order of the hypotheses in the FrequentistCalculator: FrequentistCalculator(*data, *sbModel, *bModel); //ALT=S+B, NULL=B-only
After having looked carefully at what the ProfileLikelihoodTestStat class is doing, I don’t think that only doing the thing above would be entirely correct. For discovery studies we want to evaluate the so-called q_0 test statistic, see for instance (12) in G.Cowan et.al. Hence, no matter if the pseudo-data to construct the test statistic distributions has been generated using \mu=0 (B-only) or \mu=\mu_test (S+B), the evaluation of the test statistic should always be done considering \mu=0. Correct me if I’m wrong, but I think that right now this method always evaluates the test statistic assuming \mu=\mu_test (for the conditional likelihood).
I think that this functionality could be simply added to the ProfileLikelihoodTestStat class by adding something like the following lines of code close to the beginning of ProfileLikelihoodTestStat::EvaluateProfileLikelihood:
In this way, the value for the parameter of interest used in the evaluation of the conditional likelihood will always be 0.
Am I missing anything obvious in here? This is the only way I have found to evaluate the test statistic for discovery q_0, but it requires modifying one of the original RooStats classes.
The tested value of \mu to compute the profiled likelihood ratio is provided by the used by defining a snapshot in the ModelConfig with the desired value (e.g. \mu = 0).
I have already provided a ModelConfig with a snapshop in which the POI is 0, and I think that it has been correctly picked up by the code. My concern is that I am not entirely sure that the ProfileLikelihoodTestStat class is actually using \mu=0 for the evaluation of the likelihood (as opposed to pseudo-data generation).
From lines 104 and 164 of ProfileLikelihoodTestStat.cxx I get the impression that the EvaluateProfileLikelihood function is taking the initial value of the POI to evaluate the conditional likelihood. Printing the variable initial_mu_value on screen tells me that in the evaluation of the test statistic under both the NULL (pseudo-data(\mu=0)) and ALT (pseudo-data(\mu=\mu_test)) models, the tested \mu is always equal to “\mu_test”, rather than 0.
I have also tried to use the B-only model in the definition of the test statictic object, like in the example you sent me: ProfileLikelihoodTestStatMod profll(*bModel->GetPdf())
However, this did not change the value of initial_mu_value.
Any thoughts of what might be happening? Also, could you confirm that for both the NULL and ALT distributions you are really testing \mu=0, and not \mu=\mu_test, for each step in the inverted hypothesis testing?
It is correct that the ProfileLikelihoodTestStat https://root.cern.ch/doc/v606/ProfileLikelihoodTestStat_8cxx_source.html class will use the provided test value of mu to evaluate the test statistic. But when computing a discovery significance, this test value is the one of the B model (i.e. mu = 0).
In case of a discovery significance there is no scan to different value of mu, only one test (Null hypothesis mu = 0 vs Alt hypothesis mu = standard value).
I am sure the test statistic , shown for example in the plot in
is evaluated at mu = 0 for the B toys otherwise you cannot get a half central chi2 dsitribution for the red histogram. If it is evaluated at a different test value than the one used to generate toys you will get a non-central chi-square, as it is the case for the blue histogram (S+B).