2
$\begingroup$

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]

scatter plot with tilted ellipse 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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.