[R] points() rendering points outside of input

Adam Hyland protonk at gmail.com
Fri Mar 18 21:28:44 CET 2011


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: protonk at gmail.com
web: http://en.wikipedia.org/wiki/User:Protonk



More information about the R-help mailing list