public
Last active

More fooling around with pi

  • Download Gist
pi day plus 1.R
R
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);

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.