Re: [R] points() rendering points outside of input

From: Peter Ehlers <ehlers_at_ucalgary.ca>
Date: Fri, 18 Mar 2011 17:02:46 -0700

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.



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. Received on Sat 19 Mar 2011 - 00:05:28 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sat 19 Mar 2011 - 00:50:23 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive