Home | News | Documentation | Download

Using TVirtualFitter to change error definition "UP"


#1

Hello,

I am fitting a 1D histogram with a TF1, using the TH1 method “Fit”. In the fit options, I specify a likelihood fit with Minos error estimation. What I would like is the 90% confidence limit for the returned parameter (as opposed to the 1-sigma error that is returned by default).

I am steering the fit using the TVirtualFitter class (ie. I set the default minimizer, tolerance, etc this way). Here, I use TVirtualFitter::SetErrorDef(double err) to try and change the value my minimizer uses when supplying the error in the fit parameter. However, what I find is that regardless of the value of err I supply, the returned error of my parameter never changes (I obtain it using the TVirtualFitter method GetErrors).

Looking at a scan of the likelihood space, I should certainly be seeing a different error for err=1 versus err=3, say.

Either I am doing something wrong, or misunderstand what is being represented by GetErrors (I assume this method returns the upper and lower errors supplied by the Minos routine, and that these values should be affected by the value I set err to be).

A simple version of my code follows that shows the problem. Should I be doing this in a different way completely (keeping in mind the problem I am solving is more complicated)? Any help would be great.

Thanks,

Logan

double func(double *x, double *par){
    double f = 0.;
    f = par[0]*pow(x[0],2);
    return f;
}

void test_fit(double errdef = 0.5){

TVirtualFitter::SetDefaultFitter("Minuit2");
TVirtualFitter::SetMaxIterations(7000);
TVirtualFitter::SetPrecision(1.e-8);
TVirtualFitter::SetErrorDef(errdef);
double d[21] = {102.1,79.3,64.9,50.6,37.2,23.8,13.5,9.0,4.1,0.8,0.1,1.2,3.8,9.9,
15.3,24.9,35.7,50.2,63.8,82.5,99.6};
TH1D *data = new TH1D("data","data",21,-10.,10.);
for(int i=1; i<=data->GetNbinsX(); i++){
   data->SetBinContent(i,d[i-1]);
}
data->Sumw2();

TF1 *fdata = new TF1("fdata",func,-10.,10.,1);
fdata->SetParameter(0,1.);

TFitResultPtr ptr = data->Fit(fdata,"WLVES","",-10.,10.);
int status = ptr;
if(status != 0) return;

TVirtualFitter *fit = TVirtualFitter::GetFitter();

double eup, elow, epar, globcc;
fit->GetErrors(0,eup,elow,epar,globcc);
double p = fit->GetParameter(0);

TGraph *g = new TGraph(1000);
TBackCompFitter *bfit = (TBackCompFitter*)fit;
bfit->Scan(0,g,p+10.*elow,p+10.*eup);

TCanvas *c1 = new TCanvas("c1","c1");
TCanvas *c2 = new TCanvas("c2","c2");

c1->cd(); data->Draw();
c2->cd(); g->Draw("a*");

return;
}

#2

Hi,

The problem is that you are using a weighted likelihood fit, option “WL”. In that case one can compute the error only at 1-sigma level using the inverse of the covariance matrix.
By the way, why are you using a likelihood fit ? Your data points are counts (or weighted counts) ?
If this is the case you should also set the errors as the square root of the sum of the weight square for each bin.

Lorenzo


#3

Hi Lorenzo,

Thanks for the reply. Is the answer simply that I cannot change the error definition if doing a likelihood fit (even unweighted)?

Lets say I create a histogram filled with random values pulled from a Gaussian, where the bin contents are counts. After filling the histogram, I do Sumw2(), which I believe then sets the bin errors to the square root of the bin contents (because I have not filled the histogram using weights). Now, I want to perform a likelihood fit of a Gaussian function to the histogram, and would like the 90% confidence limit, say, on the fit parameters, using h->Fit(func,“LE”), having set the error definition to ~2.3 (because it is a two parameter likelihood fit). Because I ask the fit to use the Minos error method, I assume it does not rely completely on the inverse of the covariance matrix (from the MINUIT manual), so it can find asymmetric errors. What I find is that only the 1-sigma errors are returned. If I do the default chisquared fit, then the 90% errors are returned. I’m not sure why this is the case.

In my real problem, I am fitting a binned pdf to toy MC data that is scaled to an expected rate, so my bins are not necessarily integer height (although, this could be changed). Also, it could be that my bins have few events, thus why I am doing a likelihood fit over the chisquared fit.

Again, I appreciate your help!

Logan


#4

Hi Logan,

[quote]
I cannot change the error definition if doing a likelihood fit (even unweighted)?[/quote]

I should have been more precise in my last post. You can change the error definition in a likelihood fit, if is weighted or unweighted, but at the moment (due to a bug) the errors in the weighted case are not scaled.
One thing you cannot do also right now, in a weighted likelihood fit only, is compute the Minos error, thus find the error from the likelihood function values.
You can however compute them from the covariance matrix and scale them to the value you want.
I will fix in the code now so the scaling of this error is done automatically.

If you have an histogram filled with counts and scaled with some weights you can use the weighted likelihood method. But you need to call Sumw2() and also not just SetBinContent() but also SetBinError() on your histogram, since the errors are not anymore the square root of the content.
Note that if you call, as in your code, only Sumw2() and SetBinContent() the errors are set to zero.

Best Regards

Lorenzo


#5

Ahh… Thanks for clarifying that! And thanks for the additional suggestions. I look forward to the update!

Logan


#6

This problem should be fixed now in both the trunk and 5.34 patches.
Thank you Logan for reporting it

Lorenzo


#7

[quote]Hi Logan,

[quote]
I cannot change the error definition if doing a likelihood fit (even unweighted)?[/quote]

I should have been more precise in my last post. You can change the error definition in a likelihood fit, if is weighted or unweighted, but at the moment (due to a bug) the errors in the weighted case are not scaled.
One thing you cannot do also right now, in a weighted likelihood fit only, is compute the Minos error, thus find the error from the likelihood function values.
You can however compute them from the covariance matrix and scale them to the value you want.
I will fix in the code now so the scaling of this error is done automatically.

If you have an histogram filled with counts and scaled with some weights you can use the weighted likelihood method. But you need to call Sumw2() and also not just SetBinContent() but also SetBinError() on your histogram, since the errors are not anymore the square root of the content.
Note that if you call, as in your code, only Sumw2() and SetBinContent() the errors are set to zero.

Best Regards

Lorenzo[/quote]

Hey moneta, I’m starting to dig into this now, you made a bit of sense to me but I’m a newbie so I’m trying to find some kind of help manual