15 Things to do

The code is currently only able to handle diallelic loci. With multiallelic loci, say $k$ alleles, one needs the covariances between frequencies of different alleles, between different subpopulations. These are straightforward to derive but the resulting covariance matrices involved would becomes much bigger, and would differ between different loci making the code rather cumbersome to implement. I also haven't figured out how to simulate the conditional distribution conditioned, not on a single allele frequency mean, but on a vector of means $(\bar p_1, \bar p_2, \dots, \bar p_{k-1})$.