Drawing convergence ellipse in 2d

[Georg_T]

Dear Rooters,
I would like to draw a convergence ellipse for a certain intervall (for example 1,2,3 sigma) in 2 dimensions on a set of x/y-coordinates. Let me provide you most important part of the code:

int nevt = tree->Draw("y:x","","goff");
double *Xt = tree->GetV2();
double *Yt = tree->GetV1();

double_t x0 = TMath::Mean(nevt,Xt);
double_t sx = TMath::StdDev(nevt,Xt);
double_t y0 = TMath::Mean(nevt,Yt);
double_t sy = TMath::StdDev(nevt,Yt);

TF2 *f = new TF2("f", "ROOT::Math::bigaussian_pdf(x, y, [1], [3], [4], [0], [2])", -0.1, 0.1, -0.1, 0.1);
f->SetParNames("x0", "sigmax", "y0", "sigmay", "rho"); // note: -1 < rho < 1
f->SetParameter(0, x0);
f->SetParameter(1, sx);
f->SetParameter(2, y0);
f->SetParameter(3, sy);
f->SetParameter(4, 0.01);
f->SetParLimits(4,-1,1);

//[Doing the fit routine... not important here... I used a unbinned likelihood fit here... but this is not important....this part is working]

f->SetNpx(100); f->SetNpy(100);
x0 = f->GetParameter(0);
y0 = f->GetParameter(2);
Double_t sigmax = f->GetParameter(1);
Double_t sigmay = f->GetParameter(3);
Double_t rho = f->GetParameter(4);
double factor = 1./(2*3.1415*sigmax*sigmay*sqrt(1-(rho*rho)));

//Here is the point... which parameters do I have to use for a certain intervall?
double contours[] = {0.0001*factor, 0.001*factor, 0.01*factor, 0.1*factor,1*factor};
f->SetContour(5, contours);

f->Draw("cont1z");
tree->Draw("y:x","","same *");

As you can see, I can draw the ellipses from the contour-plot, but I’m not sure how to calculate a certain z-height for my wished confidence intervall…

Of course I can setup the convergence matrix, calculate the eigenvectors and so on, but there is no easy way to draw the ellipse either…
any hint?

Additional info:
Of course I know, that bigaussian_pdf is a pdf-function, that means that the integral inside by ellipse must be equal to the given confidence interval (e.g. 0.95)
The area of an ellipse is
A=PI*a*b.
But this doesn’t help me either, as the parameters sigmax and sigmay of the bigaussian_pdf are with respect to the coordinate system, not with respect to rho

Is there a way to calculate the principal axes from sigmax, sigmay and rho?

Thanks
Georg

_ROOT Version: 6.28/06
_Platform: Windows

[mdessole]

Hello @Georg_T,

I believe that the extrema for your interval are set according to the definition of f, in your case the interval should be [-0.1, 0.1] x [-0.1, 0.1]. changing these parameters you can change the visualization frame as well.

Cheers,
Monica

[Georg_T]

Hi @mdessole,
thank you very much for your reply.

Maybe asked the question in a confusing way. Its not about the extrema.
I would like to draw the ellipse which is covering a certain interval based on a fit (e.g. 95%)

see examples for a convergence ellipse over here:
https://cran.r-project.org/web/packages/ConfidenceEllipse/vignettes/confidence-ellipse.html

I guess the correct way would be to calculate the principal axes from the covariance matrix and then set a certain factor for the axis in order to draw the ellipse.

Unluckily, I’m not 100% sure how to do the calculation based on the matrix and I’m not sure for to calculate the factor either. This is for more a mathematical, than a root problem
Georg

[mdessole]

Hello @Georg_T,

thanks for clarifying your question. What you could do is computing the confidence level of your model (fitting) function after having it fitted to an histogram. You can do that using the method TVirtualFitter::GetFitter())->GetConfidenceIntervals(h2, cl), where h2 is your 2-dim histogram and cl is your confidence level, e.g. 0.95.
If you have an histogram and you want to plot the contour lines, you just need to do h2->Draw(“CONT2”). You may what to have a look at this tutorial.

Hope this helps!
Monica

[Georg_T]

Hi @mdessole,

Sorry, but gives me the confidence level of the fit - that is not the same.

Imagine the situation in 1d. A normal distribution with mu=0 and sigma=1. The confidence level for p=0.95 is [-1.96sigma;1.96sigma]. In this case its easy to to use GetQuantiles()-method and write some lines.

In the 2d case the general gaussian approximation is a bigaussian_pdf. By using this pdf I would like to get a contour drawn which covers 0.95 (integral of the pdf within the contour). Due to the shape of the bigaussian_pdf the contour will be an ellipse

Georg

[mdessole]

As you said yourself, this is a math/stat question rather than a ROOT one.
However, if I well understand what you what to do, you want to know what is the value of f(x,y) for (x,y) in the boundary of the confidence region R_p (i.e. the region inside the ellipse) for a given probability p.
Since the region of confidence is given by (see wikipedia)
R_p = \{v = (x,y): (v-\mu)^T\Sigma^{-1}(v-\mu) \leq \chi^2_2(p) \},
and its boundary, aka the ellipse, is
E_p = \{v = (x,y): (v-\mu)^T\Sigma^{-1}(v-\mu) = \chi^2_2(p) \}.
You have that the value of f on v=(x,y) in E_p is
f(x,y) = \frac{\exp(-\frac{1}{2}(v-\mu)^T\Sigma^{-1}(v-\mu))}{2\pi\sqrt{\det(\Sigma)}} = \frac{\exp(-\frac{1}{2} \chi^2_2(p))}{2\pi\sqrt{\det(\Sigma)}}.
This is the value you’re looking for.

[Georg_T]

@mdessole,
Yes, you are totally right!

But how do I calculate \Sigma and \chi^2 ?

Georg

[mdessole]

A simple web search could answer your questions

You can find here how to express \Sigma in the bivariate case, and here the formula for the \chi^2_k function.

[Georg_T]

@mdessole: Thank you,
yes I have seen this before. And it seems it is the only solution.

My hope was, that there is a method or something hidden in root like TF2::GetPDFConfidenceZone(double p=0.95)

Anyhow…thanks a lot
Georg

[Wile_E_Coyote]

It seems to me you want to set “some_chi2_contour_level = chi2_minimum_from_the_fit + some_delta_chi2”.

See the equation “(40.80)” and the “Table 40.2” in: PDG → The Review of Particle Physics (2024) → Mathematical Tools → Statistics

Try (note: for “p = 0.95” and “ndf = 2”, the “delta_chi2 = 5.9914645”):
TMath::ChisquareQuantile(0.95, 2)
1. - TMath::Prob(5.9914645, 2)

The macro attached below may also be helpful.

bigaus_prob.cxx (5.1 KB)