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.

Improve this question
$\endgroup$
4
  • $\begingroup$ Can you give the values of sigmaA and sigmaB used above? $\endgroup$
    – SHuisman
    Mar 6, 2020 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$
    – Illari
    Mar 15, 2020 at 5:27
  • $\begingroup$ Have you considered using BinormalDistribution? $\endgroup$
    – kirma
    Mar 15, 2020 at 6:29
  • $\begingroup$ @kirma Could you explain what you mean by that? I don't immediately see the application here. $\endgroup$
    – Illari
    Mar 20, 2020 at 2:29

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.