Navigation Menu

Skip to content

Instantly share code, notes, and snippets.

Created June 1, 2014 04:02
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save anonymous/7e80257348c0a340b364 to your computer and use it in GitHub Desktop.
Save anonymous/7e80257348c0a340b364 to your computer and use it in GitHub Desktop.
Simulation of herd immunity
herd<-function(vaccination.fraction,arrival.rate=1/30,end.of.time=365,frequent.flyers=0.01,annoying.beep=TRUE){
N<-10000
xy<-matrix(runif(2*N),nrow=N,ncol=2)
status<-1+rbinom(N,1,vaccination.fraction)
startstat<-status
par(mar=c(1,1,3,1))
#svir
myrisk<-c(.3,.08,1,1,1,1,0)
yourrisk<-c(0,0,1,1,1,1,0)
nextstat<-c(3,3,4,5,6,7,7)
color<-c(adjustcolor(c("black","green","red","red","red","red"),alpha.f=0.5),"royalblue")
sym<-c(1,1,rep(19,6))
i<-1
while(i<end.of.time){
if(runif(1)<arrival.rate) {
if(annoying.beep) cat("\007")
idx<-sample(N,1)
status[idx]<-3
startstat[idx]<-NA
}
dists<-get.knn(xy)
prob<- myrisk[status]*rowSums((dists$nn.dist<0.009)*yourrisk[status[dists$nn.index]])
prob[prob>1]<-1
change<-rbinom(N,1,prob)
change[myrisk[status]==1]<-1
status[change==1]<-nextstat[status[change==1]]
plot(xy,col=color[status],pch=sym[status],axes=FALSE,xlab="",ylab="")
xy<-(xy+matrix(rnorm(2*N,s=0.01),nrow=N,ncol=2)) %%1
if (frequent.flyers>0){
idx<-sample(N,N*frequent.flyers)
xy[idx,]<-runif(2*length(idx))
}
mtext(i,adj=0.5)
if (any(status>2 & !is.na(startstat))){
tt<-table(status>2,startstat)
mtext(paste("Infected|Not vaccinated =",round(100*tt[2,1]/(tt[2,1]+tt[1,1])),"%",sep=""),adj=0)
mtext(paste("Infected|Vaccinated =",round(100*tt[2,2]/(tt[2,2]+tt[1,2])),"%"),adj=1)
}
i<-i+1
}
}
# ~/Downloads/ffmpeg -r 3 -i Rplot%d.png -y -vcodec libx264 -pix_fmt yuv420p herd-90.mp4
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment