More fooling around with pi
#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