Simulation of herd immunity

herd.R

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