## exercise 26.3 
## Part 1


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

for(i in 1:reps){
alleles<-c(rep("red",200), rep("blue",200), rep("green",200), rep("yellow",200), rep("turquoise",200))
alleles<-sample(alleles,1000,replace=F)
for(j in 1:generations){
breeders<-sample(alleles,50,replace=F)
alleles<-rep(breeders,20)
	number_remaining[i,j]<-length(levels(as.factor(alleles)))
}
}

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)")