On 2011-03-18 13:28, Adam Hyland wrote:
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.
Yes, it's one of those things one learns by being bitten. You need to ensure that you're passing a matrix to points(). But once in a while there will be only one FALSE in 'ind' and then sample.out$loc[!sample.out$ind,] will be a vector. Just add the drop=FALSE argument (see ?"["). When sample.out$loc[!sample.out$ind,] is a two-element vector, say c(v1,v2), you get the usual plot of two points at x = c(1,2), y = c(v1,v2). That's why you see the point at (2,v2). There's another point plotted at (1,v1), but it gets lost among the other points, I suppose. It's also easy to see why this won't happen with larger samples: the chance of a single ind=FALSE becomes smaller. For security, I would also add the drop=FALSE argument to your 'ind=TRUE' points. Peter Ehlers
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.
______________________________________________ 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.
