Finding coverage using Bivariate Ellipse

Finding coverage using Bivariate Ellipse

I am struggling with a problem on estimation. In several instances I have
shown how to calculate a Bivariate Ellipse using a vector of points
generated by Bivariate Normal distribution. The code works fine except
that the coverage (the number of times the generated p is contained in the
estiamted ellipse) I am getting seems extremely low. I should also state
that the old version of R gave significantly different results compared to
the new version. I am now using R 3.0.0. Heres the code to be able to
reproduce the problem.
library(MASS)
set.seed(1234) x1<-NULL x2<-NULL k<-1 Sigma2 <-
matrix(c(.72,.57,.57,.46),2,2) Sigma2 [,1] [,2] [1,] 0.72 0.57 [2,] 0.57
0.46 rho <- Sigma2[1,2]/sqrt(Sigma2[1,1]*Sigma2[2,2])
eta<-replicate(300,mvrnorm(k, mu=c(-1.01,-2.39), Sigma2))
p1<-exp(eta)/(1+exp(eta)) # true p's n<-60
x1<-replicate(300,rbinom(k,n,p1[1,]))
x2<-replicate(300,rbinom(k,n,p1[2,]))
rate1<-x1/60 # Estimated p's rate2<-x2/60
library(car) Loading required package: nnet ell <- dataEllipse(rate1,
rate2, levels=c(0.05, 0.95))
library(sp) within<-point.in.polygon(p1[1,], p1[2,], ell$0.95[,1],
ell$0.95[,2])
mean(within) # coverage [1] 0.26