# Finding percentage of points within tilted ellipse?

I am trying to find the number of data points that I have that are within a tilted ellipse that I have generated to represent the 68% confidence region. I generated some data using:

y[a_,b_,x_] := a + b*x
sigma = 0.1;
xi = Table[i,{i,10}];
len = Length[xi];
nSim = 10^3;
simSets = Table[Table[{xi[[i]],RandomVariate[NormalDistribution[y[1,0.5,xi[[i]]],sigma]]},{i,len}],nSim];
yDat = Table[Table[simSets[[j,i,2]],{i,len}],{j,nSim}];


Then I found the $$\chi^2$$ for every simulated data set like so:

firstDeg[a_,b_,x_] := a + b*x
chiSq[a_,b_] = Table[Sum[(firstDeg[a,b,xi[[i]]]-yDat[[j,i]])^2/sigma^2,{i,len}],{j,nSim}];
mins = Table[FindMinimum[chiSq[a,b][[i]],{{a,1},{b,0.5}}],{i,nSim}];


For the $$a$$ and $$b$$ values, I collected the values that minimized the $$\chi^2$$:

listA = Table[a/.mins[[j,2,1]],{j,nSim}];
listB = Table[b/.mins[[j,2,2]],{j,nSim}];


and then tried finding the 68% confidence region by using the fact that the semi-major and minor axis are $$\sqrt{2.3}\sigma_a,\sqrt{2.3}\sigma_b$$

EDIT: $$\sigma_a,\sigma_b$$ found via:

VarA = Total[(listA-Mean[listA])^2]/Length[listA];
VarB = Total[(listB-Mean[listB])^2]/Length[listB];
sigmaA = Sqrt@VarA;
sigmaB = Sqrt@VarB;


and the orientation of the ellipse via $$\alpha=\arctan\left(\frac{v_1(b)}{v_1(a)}\right)$$, where $$v_1$$ refers to the eigenvector corresponding to the largest eigenvalue. I calculate the region for the parameters like so:

axes[delchisq_,sig_] := sig*Sqrt[delchisq]
paramVal = Table[{a,b}/.nlm[[j]]["BestFitParameters"],{j,nSim}];
paramValM = Mean@paramVal;
angle[covariance_] := ArcTan[Eigenvectors[covariance][[1,2]]/Eigenvectors[covariance][[1,1]]]

aUp = paramValM[[1]]+axes[2.3,sigmaA]
aLow = paramValM[[1]]-axes[2.3,sigmaA]
bUp = paramValM[[2]]+axes[2.3,sigmaB]
bLow = paramValM[[2]]-axes[2.3,sigmaB]
valIn = Table[If[(paramVal[[i,1]] >= aLow \[And] paramVal[[i,1]] <= aUp)\[And](paramVal[[i,2]] >= bLow \[And] paramVal[[i,2]] <= bUp),paramVal[[i]],Unevaluated[Sequence[]]],{i,Length@paramVal}];
Print["The percentage of (a,b) points within the 68% confidence region is ",N@(Length@valIn/nSim)*100,"%."]
(* The percentage of (a,b) points within the 68% confidence region is 83.9%. *)


and then I plot (I also calculated the covariance matrix "by hand", as it were, and that is to what covMat refers):

p68 = Graphics[Rotate[Circle[paramValM,{axes[2.3,sigmaA],axes[2.3,sigmaB]}],angle[covMat]],Frame->True];
valPt = Graphics[{Directive[Red],PointSize[Medium],Point[valIn]},Frame->True];
Show[p68,valPt]


Obviously, I'm not counting the points in the ellipse correctly, but beside that, I'm not counting many points outside of the ellipse, and yet I'm getting that 83.9% of my $$(a,b)$$ points are in the region. Why would that be?

Any help in correctly 1) calculating the 68% confidence region (I'd prefer not to use ["ParameterConfidenceRegion", ConfidenceLevel -> 0.68]), and 2) finding the points within the tilted ellipse would be a great help.

• Can you give the values of sigmaA and sigmaB used above? – SHuisman Mar 6 '20 at 9:52
• @SHuisman oh sure, it's been edited in. I think it would be about what one would expect, though? Maybe I'm wrong tho – Illari Mar 15 '20 at 5:27
• Have you considered using BinormalDistribution? – kirma Mar 15 '20 at 6:29
• @kirma Could you explain what you mean by that? I don't immediately see the application here. – Illari Mar 20 '20 at 2:29