Here is a program to simulate genetic drift in R (code pasted below), specifically one **capable of handling unequal sex ratio**. Sex ratio directly affects the effective population size, and this simulation can be used to show that it is the effective population size (*N*_{E}), rather than the census population size (*N*_{A}), that determines the population genetic behaviour of a population.

Here are some demonstration results, where 10 replicate populations are simulated with *N*_{E} = 360, but very different actual population sizes, and an initial allele frequency of 0.5. The left hand simulations have unequal sex ratios, while the sex ratio is equal in the right hand simulations. Both sets of simulations behave similarly, showing **that effective population size, rather than actual population size, is the determinant of population genetic behaviour.**

#####~ Genetic Drift Simulator ~#####
#######~ Robert M. Griffin ~#######
#####~ Updated March 22nd 2016 ~#####

# Pop = Replicate populations
# Gen = Generations
# NM = Male population size
# NF = Female population size
# P = Frequency of focal allele
# histo = Optional histogram

GenDriftSim = function(Pop = Pop, Gen = Gen, NM, NF, P, histo = "yes"){
P = (2*(NM+NF))*P
NE = round((4*NM*NF)/(NM+NF),0)
SR = round(NM/NF,2)
Na = NM+NF
plot(c(0,0),type = "n", main = bquote('N'[M]~'/ N'[F]~'='~.(SR)*', N'[A]~'='~.(Na)*', N'[E]~'='~.(NE)), cex.main = 1,
xlim = c(1,Gen), ylim=c(0,1), xlab = "Generations", ylab = "Fequency of focal allele")
for (i in 1:Pop){
N = NM+NF
startA = as.vector(c(rep(1, times = P),rep(0, times = (2*N)-P)))
Population = matrix(c(
c(sample(startA, size = 2*N, replace = FALSE)),
c(rep("M", times = NM), rep("F", times = NF))),
ncol = 3)
SimResults[(Gen*i)+1-Gen, 3] <<- sum(as.numeric(Population[,1:2]))/(N*2)
for(j in 1:(Gen-1)){
Population = matrix(c(
c(sample(sample(Population[(1:NM),1:2], replace = TRUE),N, replace = TRUE)),
c(sample(sample(Population[(1+NM):N,1:2], replace = TRUE),N, replace = TRUE)),
c(rep("M", times = NM), rep("F", times = NF))), ncol = 3)
SimResults[(Gen*i)+1+j-Gen, 3] <<- sum(as.numeric(Population[,1:2]))/(N*2)
}
s = (i*Gen)-Gen+1; e = i*Gen
r = as.vector(SimResults[s:e, 3])
points(r~c(1:Gen), type = "l")
}
if(histo == "yes"){SimResults[,1] = rep(1:Pop, each = Gen)
SimResults[,2] = rep(1:Gen, times = Pop)
hist(SimResults[,3][SimResults[,2]==Gen], xlim = c(0,1), cex.main = 1, main = bquote('N'[M]~'/ N'[F]~'='~.(SR)*', N'[A]~'='~.(Na)*', N'[E]~'='~.(NE)), xlab = paste0("Frequency of focal allele after ",Gen," Generations"))
}else{}
}
Pop = 10
Gen = 25
P = 0.5
histo = "no"
par(mfrow=c(1,2))
SimResults = matrix(data = NA, ncol = 3, nrow = Gen*Pop)
GenDriftSim(Pop = Pop, Gen = Gen, NM = 100, NF = 900, P = P, histo = histo)
GenDriftSim(Pop = Pop, Gen = Gen, NM = 180, NF = 180, P = P, histo = "histo")