From: Peter Ehlers <ehlers_at_ucalgary.ca>

Date: Fri, 18 Mar 2011 17:02:46 -0700

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

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.
*