Quit[];

(* The magic ((δχ)^2) numbers *)
dx2=2 InverseGammaRegularized[(M/2),0,Erf[(n/(Sqrt[2]))]]/.M->2/.n->{1,2,3};

(* Assume a covariance matrix, make a (χ^2) and make some mock data *)
Cij={{1,1/10},{1/10,1}};
chi2={x,y}.Inverse[Cij].{x,y};
data=RandomVariate[MultinormalDistribution[{0,0},Cij],{100000}];

(* Find the empirical distribution of the data, its PDF and Quantiles at 1,2,3σ *)
skd=SmoothKernelDistribution[data];
pdf=PDF[skd,data];
cons=Quantile[pdf,1-Erf[(nn/(Sqrt[2]))]/.nn->{1,2,3}//N];
cenp=data//Mean;

(* Plot the contours for the empirical and analytic distribution using the (χ^2) to see if they match *)
op=0.5;
cols1={Transparent,{Opacity[op],Hue[0.6,0.8,1]},{Opacity[op],Hue[0.55,0.8,1]},{Opacity[op],Hue[0.5,0.8,1]}};
cols2={Transparent,{Opacity[op],Hue[0.1,0.8,1]},{Opacity[op],Hue[0.15,0.8,1]},{Opacity[op],Hue[0.2,0.8,1]}}//Reverse;
con1=ContourPlot[PDF[skd,{a,b}],{a, -4,4},{b, -4, 4},Contours->cons,ContourLines->True,ContourShading->cols1,ContourStyle->Red,FrameLabel->{"x","y"},BaseStyle->FontSize->18,PlotRange->All,FrameStyle->Black,ClippingStyle->None,Epilog->{Black,Point[cenp]},FrameStyle->Black];
con2=ContourPlot[chi2,{x, -4,4},{y, -4, 4},Contours->dx2//Evaluate,ContourLines->None,ContourShading->cols2,ContourStyle->Black,FrameLabel->{"x","y"},BaseStyle->FontSize->18,PlotRange->All,FrameStyle->Black,ClippingStyle->None,Epilog->{Black,Point[cenp]},FrameStyle->Black];
{con2,con1}

(* Perfect match! *)
Show[{con2,con1}]