Bias in HistFactory fit

Roofit and RooStats
1

Dear experts,

For a while, I have been trying to troubleshoot a bias I see in my fits with HistFactory (with ROOT 6.22.06).

For the purpose of this post, I set up toys with a very simple model of a Gaussian signal (mean = 0, sigma = 1) over a uniform background, in two channels. The model for Channel 1 has signal + bkg, and Channel 2 has only background.

The toy data is generated to be background only with a yield of 10k in each channel per toy.

The fit is done between -5 and 5, with 40 bins.

The background shape in the fit is data driven using HistFactory::Sample::AddShapeFactor(). The ShapeFactor is shared between the two channels.

The signal shape is kept same from toy to toy. The generated bkg only data and background template are varied from toy to toy (using TH1::FillRandom()).

When I run the fit over 1000 toys, I see a bias in the distribution of the fitted POI value. The numerical mean of the distribution is -3.9 +/- 3.3.

I see this bias in many instances of this fit, with different signal widths, different number of bins, different generated yields, and I am not sure of why it is there.

I would appreciate any advice on what I’m doing wrong. I have attached a script that should reproduce what I am doing.

Thanks,
Arvind.

PS: When I run this script, I also see the following warning for each toy
WARNING: Can't find parameter of interest: nSig in Workspace. Not setting in ModelConfig.

I think this come from the fact that my POI (nSig) is in only one channel and not the other, so HistFactory complains a bit.

troubleShootHF_2Channel.C (8.8 KB)

2

Dear experts,

Sorry for bothering you again, but have one of you had a chance to look at this?

Since last posting, I have also seen bias in the following variations of the above:

  1. Single channel fit with a known flat background (i.e. no ShapeFactor)
  2. Using TRandom3 instead of FillRandom to generate the background data and templates.
  3. Using RooFit with HistPdfs instead of HistFactory to do the fit.

I can also say that the fitted values I get from HistFactory are very consistent with those that I get from RooFit+HistPdf.

So if two different fitting methods both show bias, this leads me to think that something must be wrong with the way I am generating the bkg only data and templates, but I can’t figure out what.

Any assistance would be much appreciated,
Arvind.

3

Hi,

Sorry for the late reply. I will look at it later today
Cheers

Lorenzo

1 Like
4

Dear @avenkate,

I see one problem with how you generate the Gaussian signal template, you are filling the template with TH1::FillRandom, so statistical fluctuations in the template creation can cause a bias.

You could create a function to create the template such that it represents the expected number of counts in each bin (Asimov data). You can write a function that does this by copy-pasting TH1::FillRandom and change the line with the random number generator.

void TH1FillAsimov(TH1 & th1, const char *fname, Int_t ntimes)
{
   Int_t bin, binx, ibin, loop;
   Double_t r1, x;
   //   - Search for fname in the list of ROOT defined functions
   TF1 *f1 = (TF1*)gROOT->GetFunction(fname);
   if (!f1) { Error("FillRandom", "Unknown function: %s",fname); return; }

   //   - Allocate temporary space to store the integral and compute integral
    
   //TAxis * xAxis = &fXaxis;
   TAxis * xAxis = th1.GetXaxis();

   // in case axis of histogram is not defined use the function axis
   if (xAxis->GetXmax() <= xAxis->GetXmin()) {
      Double_t xmin,xmax;
      f1->GetRange(xmin,xmax);
      Info("FillRandom","Using function axis and range [%g,%g]",xmin, xmax);
      xAxis = f1->GetHistogram()->GetXaxis();
   }                                

   Int_t first  = xAxis->GetFirst();
   Int_t last   = xAxis->GetLast();
   Int_t nbinsx = last-first+1;

   std::vector<Double_t> integral(nbinsx+1);
   integral[0] = 0;
   for (binx=1;binx<=nbinsx;binx++) {
      Double_t fint = f1->Integral(xAxis->GetBinLowEdge(binx+first-1),xAxis->GetBinUpEdge(binx+first-1), 0.);
      integral[binx] = integral[binx-1] + fint;
   }

   //   - Normalize integral to 1
   if (integral[nbinsx] == 0 ) {
      Error("FillRandom", "Integral = zero"); return;
   }
   auto normInv = 1./integral[nbinsx];
   for (bin=1;bin<=nbinsx;bin++)  integral[bin] *= normInv;

   //   --------------Start main loop ntimes
   for (loop=0;loop<ntimes;loop++) {
      r1 = static_cast<double>(loop) / ntimes; // this is different compared to TH1::FillRandom
      ibin = TMath::BinarySearch(nbinsx,&integral[0],r1);
      //binx = 1 + ibin;
      //x    = xAxis->GetBinCenter(binx); //this is not OK when SetBuffer is used
      x    = xAxis->GetBinLowEdge(ibin+first)
             +xAxis->GetBinWidth(ibin+first)*(r1-integral[ibin])/(integral[ibin+1] - integral[ibin]);
      th1.Fill(x);
   }
}

For the background, you already create such an Asimov template with your makeFlatHisto.

With an Asimov template also for the signal, I don’t see any bias for nSig in 1000 fitting rounds. The mean of nSig is -4.06 +/- 3.34, and more importantly: the mean of the pull distribution nSig/nSig_err is -0.028 +/- 0.0319. The sigma of the pull distribution is 1.0079. So there is no bias in your HistFactory fit.

In fact, the mean -3.9 +/- 3.3 you obtained is also compatible with zero, and to be more correct you should look at the mean of the pull distribution anyway. Hence I think that the fluctuations in your signal template are not a big deal. But to be correct I think it’s still important to use the Asimov templates.

I hope this clarifies what you saw, and if not please feel free to ask again!

Cheers,
Jonas

1 Like
5

Hi Jonas,

Thanks so much for getting back to me. I am still processing your reply. But I wanted to clarify some things you said:

  1. As long as the mean and sigma of the pull distribution are consistent with 0 and 1 respectively, then we can conclude that there is no significant bias, even if the mean of the POI is (say more than 1 sigma) away from zero.

  2. Further, in each toy, the signal and background templates should represent exactly the expected distributions. There should be no fluctuations in these templates from toy to toy.

Have I understood these two things correctly?

Cheers,
Arvind

6

Hi Arvind, yes you summarized that correctly! You can also think about your perceived “bias” from another angle. In your script, you have set genYield = 10000, which is used both as signal and background yield in channel 1.

Since you don’t include signal in your data (null hypethesis), your fitted nSig will be roughly of the order of magnitude of the background fluctuations, which is sqrt(10000) = 100. So if your POI is expected to be of the order of 100 away from zero, you should not be surprised that it had a mean of roughly 5 in our toy experiments.

1 Like
7

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