As a followup to pi-day, I attempted to make a .gif of a simulation
based estimation of pi by plotting points inside a single quadrant of
a circle (a la http://www.drewconway.com/zia/?p=2667 ). When
rendering the individual x,y pairs with points() I intermittently see
points crop up around (2,0.5) but the input values for x and y are
bounded between 0 and 1.
square<-structure(c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0), .Dim = c(5L, 2L));
library(animation)
base.plot<- function() {
plot(-2:2,-2:2,type="n",asp=1)
lines(square)
symbols (0, 0, circles=1, inches=FALSE, add=TRUE)
}
pi.ratio<- function(n) {
x<- runif(n, min=0,max=1)
y<- runif(n, min=0,max=1)
pi.pts<- cbind(x,y)
return(list(est=4*sum(x*x + y*y <= 1.0)/n, ind=x*x + y*y <= 1.0,
loc=pi.pts))
}
pi.est<- function(iter,n) {
sum.pi<- stor.rat<- numeric(0)
for (i in 1:iter) {
sample.out<- pi.ratio(n)
stor.rat[i]<- sample.out$est
sum.pi[i]<- sum(stor.rat[1:i])/i
base.plot()
points(sample.out$loc[sample.out$ind,],col=2,pch=20,cex=0.8)
points(sample.out$loc[!sample.out$ind,],col=4,pch=20,cex=0.8)
}
}
saveMovie(pi.est(50,20),interval = 0.03, movie.name =
"pianim.gif",outdir = getwd(), width = 600, height = 600);
If you don't have animation installed and don't need to see a .gif you
can just pull base.plot() out of pi.est() and render successive loops
on top of each other. After a few dozen iterations you will see a
point plotted well outside of the constraints for x,y. I'm not sure
what causes this behavior. Running pi.ratio() for very large sample
sizes never returns an x value greater than or equal to one (which
accords with the documentation for runif()).
I'm pretty sure this is some subtle (or not so subtle) problem
stemming from my inexperience and not a bug. :) I would appreciate
any help that can be offered.
Thanks
--
Adam Connors Hyland
email: prot...@gmail.com
web: http://en.wikipedia.org/wiki/User:Protonk
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
