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.