Created
March 15, 2011 20:58
-
-
Save Protonk/871478 to your computer and use it in GitHub Desktop.
More fooling around with pi
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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); |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment