[R] points() rendering points outside of input

[Peter Ehlers]

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