set.seed(123)
xxe <- runif(10000,2,10)
yye <- runif(10000,3,9)

radioa <- 4
radiob <- 3
c <- sqrt(radioa^2-radiob^2)
plot(xxe,yye, cex=0.3, xlim = c(2,10), ylim = c(3,9), yaxt = "n")
axis(2, at = c(3, 6, 9))
points(6,6, cex=2, col="blue", pch=16)
points(6-radioa,6, col="red");points(radioa+6,6, col="red")
points(6,6-radiob, col="red");points(6, 6+radiob, col="red")
points(6+c,6, col="red");points(6-c,6, col="red")

foco1<-c(6+c,6)
foco2<-c(6-c,6)

distancia <- function (x1,y1,foco1,foco2){
  d1=sqrt((x1-foco1[1])^2+(y1-foco1[2])^2)
  d2=sqrt((x1-foco2[1])^2+(y1-foco2[2])^2)
  d1+d2
}

count_e=0
for (i in 1:10000){
  if (distancia(xxe[i],yye[i],foco1,foco2) < 2*radioa){
    count_e=count_e + 1
    points(xxe[i], yye[i], col= "green", pch=16, cex=0.3)
  }
}

count_e
## [1] 7894
puntos_e <- count_e/length(xxe)
area_rectangulo <- 8*6
area_estimada_e <- puntos_e*area_rectangulo
pi_estimado_e <- area_estimada_e/(radioa*radiob)
data.frame(area_estimada_e, pi_estimado_e)
##   area_estimada_e pi_estimado_e
## 1         37.8912        3.1576