Protonk / pi day plus 1.R

Created March 15, 2011


    
    View pi day plus 1.R 
  
  
      pi day plus 1.R
    
    
        Raw
      
    
    File suppressed. Click to show.
  
  
            1
            2
            3
            4
            5
            6
            7
            8
            9
            10
            11
            12
            13
            14
            15
            16
            17
            18
            19
            20
            21
            22
            23
            24
            25
            26
            27
            28
            29
            30
            31
            32
            33
            34
            35
            36
            37
            38
            39
            40
            41
            42
            43
            44
            45
            46
            47
            48
            49
            50
            51
            52
            53
            54
            55
            56
            57
            58
            59
            60
            61
            62
            63
            64
            65
            66
            67
            68
            69
            70
            71
            72
            73
            74
            75
            76
            77
            78
            79
            80
          
          
            #Same idea as yesterday, but in one function with two arguments
pi.est<- function(iter,n) {
	pi.ratio<- function(n) {
		x<- runif(n, min=-1,max=1)
		y<- runif(n, min=-1,max=1)
		return(4*sum(x^2 + y^2 <= 1.0)/n)
		}
	sum.pi<-stor.rat<-numeric(0)
	for (i in 1:iter) {
		stor.rat[i]<- pi.ratio(n)
		sum.pi[i]<- sum(stor.rat[1:i])/i
		}
	return(sum.pi)
	}
 
#Takes the function pi.est() and runs it for 10,100, and 1000 within group random
#variables so we can compare convergence
mapply.test<-mapply(pi.est,2000,c(10,100,1000))
colnames(mapply.test)<-c("10","100","1000")
convergence<-as.data.frame(mapply.test)
convergence$iters<- 1:2000
library(ggplot2)
gg.converge<-melt(convergence,id.vars="iters")
colnames(gg.converge)<-c("iters","RVs","Estimate")
conv.plot<-ggplot(data=gg.converge,aes(x=iters,y=Estimate))
conv.plot + geom_line(aes(colour=RVs)) + opts(title="Estimate of pi with different \n within group sample sizes")
 
#How much does the estimate change over iterations?
conv.test<-diff(abs(mapply.test - pi))
diff.conv<-as.data.frame(conv.test)
diff.conv$iters<- 1:1999
gg.diff<- melt(diff.conv,id.vars="iters")
colnames(gg.diff)<- c("iters","RVs","Diff")
diff.plot<- ggplot(data=gg.diff,aes(x=iters,y=Diff))
diff.plot + geom_line(aes(colour=RVs))+  opts(title="Rate of convergence with different \n within group sample sizes")
 
####Animated breakdown of the process above #####
 
#Nearly empty plot.  The only interesting bit is the tacked on y-axis for the running average sparkline
base.plot<- function() {
	plot(-2:2, -2:2, type="n", asp=1, ann=FALSE, axes='FALSE', frame.plot=TRUE)
	title(main="A peek inside the \n estimation of pi")
	axis(side=2,at=c(seq(0.25,2,by=0.25)),labels=c(as.character(seq(0.25,2,by=0.25)+2)),par(cex=0.8,las=1))
	lines(structure(c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0), .Dim = c(5L, 2L)))
	symbols (0, 0, circles=1, inches=FALSE, add=TRUE) 
	}
 
#Same as the function above with 0,1 as bounds.  I return three objects because I want to get at the 
#x,y coords, the location and the ratio all at once (outside the for() loop).
 
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))
	}
	
#Nothing too dramatic, though the ",drop=FALSE" argument prevents the red and blue matrices from 
#accidentally becoming vectors.
 
pi.est<- function(iter,n) {
	sum.pi<- stor.rat<- numeric(0)
	x.len<-seq(-2,2,length.out=iter)
	for (i in 1:iter) {
		sample.out<- pi.ratio(n)
		red<- as.matrix(sample.out$loc[sample.out$ind, , drop=FALSE])
		blue<- as.matrix(sample.out$loc[!sample.out$ind, , drop=FALSE])
		stor.rat[i]<- sample.out$est
		sum.pi[i]<- sum(stor.rat[1:i])/i
		base.plot()
		points(red,col=2,pch=20,cex=0.8)
		points(blue,col=4,pch=20,cex=0.8)
		text(-1.5,0.8,labels=sprintf("%f",sum.pi[i]),col=3)
		lines(x.len[1:i],sum.pi[1:i]-2,col=3)
		Sys.sleep(0.02);
		}
	}
 
library(animation)
saveMovie(pi.est(50,20),movie.name = "pianim.gif",outdir = getwd(), width = 600, height = 600);