Minuit failed in fitting a 2d histogram of correlated variables. Thus I projected along x and y, and from the fit of the projections I obtained sigmax and sigmay. Now I would like to plot 1 sigma and 2 sigma contours. Is there any simple way to do that? I know that Minuit provides such a functionality, but as far as I know it works for a 2d fit.
If you want the correct contour you need to make a 2D fit. Otherwise your contour is just a simple ellipse (not rotated) with axis sigma and sigma.
Thank you for answering.
Yes, I know that this is not the correct ellipse, but I thought that there was a way to get the correct ellipse through a transformation involving the covariance matrix.
I would prefer to go through a 2d fit, but this does not work. I hardly understand why since the projections are perfect gaussians.
Yes you can probably use the covariance matrix estimated from just the data points and build a 2d-gaussian and draw the contour of that function.
You can do this in ROOT by making a TF2 object and then set the correct contour levels using TF2.SetContour.
However, it is strange that the fit is not working. What are you fitting exactly the data with ? A 2d-gaussian ?
Actually I wonder if I have not already filled the 1 sigma contour.
I may be worth explaining the context:
I am doing n experiments to measure mass and width of the W from the measurement of the cross-section at threshold. To a pair (mass,width) corresponds a cross-section. Then I compare the pseudo-measured cross-section to the cross-sections in a table, and I build a 2d-chi2. For each experiment, I fit the projections of this 2d-chi2. Then I fill a 2d-histogram with (mass, width) for each experiment, such that I have n entries in the histogram. The result is an ellipse, showing the correlation of mass and width, the mean (x,y), i.e. the center of the ellipse has the coordinates that have been obtained in the 1d-fits, and the standard deviations are also those from the 1d-fits.
Then I tried to fit interactively with the fit panel, using a 2d-gaussian. This is probably stupid, as it fails.
In those conditions, if I get you right, TF2.SetContour would do the job. Did I get you right?
I had not used SetContour before. I am not sure about how to set the contour levels, since, unless I am wrong, a given contour will enclose the region where the bins are below/above e a given threshold (i.e., if I set contour levewl to 10, I will get all the bins with a content<10), not the region containing a given percentage of the population (this is what I would need actually).
As the 2D fit is not supposed to fail, I tried bigaus that is supposed to be a member of TFormula (according to what I understood from the forum). But when trying to use it, I get the message “no member called bigaus in TFormula”. Anybody having any idea about what is wrong?
There must be a simple way to fit and get the nsigma contours of a simple bivariate normal distribution. Any good advice to get rid of this?
Thanks for your help.
Try “xygaus” (parameters are “const”, “meanx”, “sigmax”, “meany”, “sigmay”):
TF2 *f = new TF2("f", "xygaus", -5, 6, -7, 8);
f->SetParameters(10, 1, 3, -2, 2);
f->Draw("colz"); or try this “bivariate-gaussian” (parameters are “const”, “meanx”, “sigmax”, “meany”, “sigmay”, “correlation”):
TF2 *f = new TF2("f", "*TMath::Exp(-0.5/(1.0-*)*(TMath::Power((x-)/,2)+TMath::Power((y-)/,2)-2.0**(x-)/*(y-)/))", -5, 6, -7, 8);
f->SetParameters(10, 1, 3, -2, 2, 0.5);
P.S. I assume “bigaus” exists only in ROOT master (as of today).
Thank you for the explanation.
For the fitting of the histogram, which function did you use from the fitpanel ? You should use the “bigaus” function, available in the latest ROOT versions and you might need to set some sensible initial parameters before fitting.
I attached here an example of fitting a 2d histogram. But, yes bigaus exists only from ROOT 6.07.06 and the master.
exampleBigaus.C (896 Bytes)
Dear Pepe and Lorenzo,
Many thanks! Before your last post, Lorenzo, I tried the suggestion of Pepe and now IT WORKS! The fit is excellent and now I have to draw the 1 and 2 sigma contours, which should be easy then, I guess.
Thanks again for your kind and efficient help,
BTW. As Lorenzo reminds … for your fitting function, your should “set some sensible initial parameters before fitting”, otherwise the fit may fail (especially if you “manually” create your own “bivariate-gaussian” function).
I forgot to tell you that my Root version is 6.06/02
No problem. Both solutions, mine (“xygaus” and “bivariate-gaussian”) and Lorenzo’s (“exampleBigaus.C”), will work in any ROOT version.
It looks like my message before that giving my ROOT version was not taken.
I was saying that the fit is now perfect, but I get a strange message when trying to draw the 1 sigma and 2 sigma contours: error: use of undeclared identifier 'gMinuit’
Did I miss anything trivial? I checked with an example given by René on the forum, and I can’t find what’s wrong.
Here is the short piece of code:
TF2* fitEllipse = new TF2("fitEllipse","*TMath::Exp(-0.5/(1.0-*)*(TMath::Power((x-)/,2)+TMath::Power((y-)/,2)-2.0**(x-)/*(y-)/))",80.383,80.388,2.047,2.053); fitEllipse->SetParameters(1., 80.3853,0.0004891,2.0495,0.00062663,-0.698657); hMassVsWidth->Fit(fitEllipse); gMinuit->SetErrorDef(1); TGraph *gr1 = (TGraph*)gMinuit->Contour(80,1,3); gr1->SetFillColor(42); gr1->Draw("alf"); gMinuit->SetErrorDef(4); TGraph *gr2 = (TGraph*)gMinuit->Contour(80,1,3); gr2->SetFillColor(38); gr1->Draw("lf");
Attached is an example making the contour.
exampleBigaus.C (1.43 KB)
Thank you for this piece of code.
I tried it, but got this error:
error: use of undeclared identifier 'TBackCompFitter
It looks like I am missing an include statement.
Here are my include statements:
If you are compiling the macro, yes, you might need to add
Thank you! This works formally.
However I need to show on the plot not sigmax vs sigmay, but mass vs width (as in the fitted plot attached), thus I guess I have to replace
par1=2 by par1=1, and par2=4 by par2=3.
Doing that I obtain the nice contour plot attached, but clearly the contours are not consistent with the result of the fit:
FCN=11215.4 FROM MIGRAD STATUS=CONVERGED 301 CALLS 302 TOTAL
EDM=1.96703e-06 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 p0 2.16972e+00 1.95199e-02 5.69734e-04 -2.84023e-04
2 p1 8.03853e+01 8.43708e-06 3.83307e-05 1.98930e+02
3 p2 9.07349e-04 1.39004e-05 4.48770e-07 2.34219e+01
4 p3 2.04995e+00 1.04738e-05 9.77493e-07 -3.15825e+01
5 p4 1.13604e-03 1.73521e-05 5.64954e-07 -8.04526e+00
6 p5 -7.03939e-01 9.67108e-03 3.54909e-04 -1.87776e-02
Why do you say it is too small ? For me it looks consistent with the parabolic error you get from the fit.
Maybe you don’t want the contour on the fitted values of M and width (which is the basically the average over all your events), but just the spread of the M and Width. In this case you need to draw the contours of the fitted 2d gaussian functions not from the covariance matrix from the fit
Yes, I probably did not express what I wanted correctly.
What I want to draw is the expected mass and width, and the expected statistical errors from the simultaneous fit, i.e. basically the area laying between (mean mass +/- n sigma mass) vs (mean width +/- nsigma width) with n=1 or 2.
Then I may have to draw the fitting function with the fitted mass and width with the axes equal to 1 or 2 sigma in both directions. Is this correct?