nextgeneration <- function(M,G) { ## M[i,j] is allele (0 or 1) of individual j in generation i ## G[i,j] is the mother of individual j of generation i i <- dim(M)[[1]] J <- dim(M)[[2]] M <- rbind(M,rep(-1,J)) G <- rbind(G,rep(-1,J)) z <- 0 ## counter of individuals of genereatoin j+1 who already know their ## parent for(j in 1:J) { ## let j have offspring if(j==J) { k <- J-z } else { ## print(c(J,z,j)) k <- rbinom(1,size=J-z,p=1/(J-j+1)) } if(k>0) { ## print(c(dim(G),i,z,k)) G[i+1,(z+1):(z+k)] <- rep(j,k) M[i+1,(z+1):(z+k)] <- rep(M[i,j],k) } z <- z+k } return(list(M,G)) } driftsimulator <- function(popsize,generations,initial) { ## population size, number of generations to be simulated, ## and initial allele frequency M <- matrix(c(rep(1,initial),rep(0,popsize-initial)),nrow=1,ncol=popsize) G <- matrix(rep(NA,popsize),rep(0,popsize-initial),nrow=1,ncol=popsize) for(i in 2:generations) { L <- nextgeneration(M,G) M <- L[[1]] G <- L[[2]] } return(list(M,G)) } driftplot <- function(M,G,generations=dim(M)[[1]],stepwise=TRUE,nsteps=10) { g <- generations ## number of generations n <- dim(M)[[2]] ## population size plot(sum(M[1,])/n,xlim=c(0,g),xlab="Generation",ylim=c(0,5.4), ylab="Allele frequency",pch=16,yaxt="n") axis(2,c(0,1),c(0,1)) abline(h=c(0,1)) points(rep(1,n),1.3+4*(1:n-1)/n,col=M[1,]*2,pch=16) points(rep(1,n),1.3+4*(1:n-1)/n) for(gg in 2:g) { if(stepwise && gg