drift<-function(N,G,M){
# N = number of genes in population
# G = number of generations of simulation
# M = number of A alleles in population
pop<-matrix(nrow=N, ncol=G)
# Holds resulting alleles, each generation
prop<-rep(NA,G)
# Holds proportion of allele A in each generation
prop[1]<-M/N
#Initialize the first generation
pop[,1]<-c(rep(1,M),rep(0,(N-M)))
for(j in 2:G){ #looping over generations
for (i in 1:N){
pop[i,j]<-sample(pop[,j-1],1,replace=TRUE)
}
##############################################
# Alternative code to replace interior loop: #
# pop[,j]<-sample(pop[,j-1],replace=TRUE) #
##############################################
prop[j]=sum(pop[,j])/N
}
return(prop)
}