Simple model of population genetics

GenPop.r

gfreq = function(p) list(AA=p^2,Aa=2*p*(1-p),aa=(1-p)^2)

mfit = function(p,fit)
{
	with(as.list(fit),{
		mf = (p^2)*wAA+(2*p*(1-p)*wAa)+((1-p)^2)*waa
		return(mf)
	})
}

pp = function(p,fit,mu)
{
	mf = mfit(p,fit)
	with(as.list(fit),{
		np = ((p^2)*wAA+p*(1-p)*(wAa)/(mf))*(1-mu)
		return(np)
	})
}

simul = function(p=0.1,s=0.1,mu=0.01,t=100)
{
	fit = list(wAA=1,wAa=1,waa=1-s)
    Results = matrix(0,ncol=6,nrow=t+1)
    Results[1,] = c(0,p,(1-p),unlist(gfreq(p)))
	for(i in 1:t)
	{
		p = pp(p,fit,mu)
		Results[(i+1),] = c(i,p,(1-p),unlist(gfreq(p)))	
	}
	Results = as.data.frame(Results)
    colnames(Results) = c('t','A','a','AA','Aa','aa')
	return(Results)
}

out = simul(t=200)
plot(out$t,out$Aa,lty=3,type='l',ylim=c(0,1),xlab='Generation',ylab='Genotype frequency')
lines(out$t,out$AA,lty=1,type='l')
lines(out$t,out$aa,lty=2,type='l')