## exercise 26.3 
## Part 1


reps<-100
generations<-500


number_remaining<-matrix(NA,nrow=reps,ncol=generations)

for(i in 1:reps){
start<-c(rep("a1",200),rep("a2",200),rep("a3",200),rep("a4",200),rep("a5",200))
a<-sample(start,1000)
allele1<-c(a[1:500])
allele2<-c(a[501:1000])
genotypes<-rbind(allele1,allele2)

for(j in 1:generations){
n<-sample(1:500,25)
temp<-sample(rep(c(genotypes[1,n],genotypes[2,n]),20),1000)
allele1<-c(temp[1:500])
allele2<-c(temp[501:1000])
genotypes<-rbind(allele1,allele2)
number_remaining[i,j]<-length(levels(as.factor(temp)))
}}

plot(NULL,ylim=c(0,5),xlim=c(1,generations),xlab="Generation",ylab="Mean number of alleles in population")

for(i in 1:generations) lines(c(i,i),c(mean(number_remaining[,i])-(1.96*sqrt(var(number_remaining[,i]))),mean(number_remaining[,i])+(1.96*sqrt(var(number_remaining[,i])))),pch=19,col="grey60")
for(i in 1:generations) points(i,mean(number_remaining[,i]),pch=16,col="blue",main="Mean alleles remaining (95% confidence intervals)")