simulation of reduction of variance and Bulmer effect

As it is so difficult to understand reduction of variance due to selection and Bulmer effect I just simulated it in julia (code below). This is a rather crude simulation,  but it gives the following plot according to expectattions (the plot using Plotly is actually quite poor).

# simulation of Bulmer effect

using LinearAlgebra
using Distributions
using Random
#using Plotly
using PyPlot

Random.seed!(1234);

nqtl=300
nanim=10000
nselect=Int(nanim/20)
h2=0.8
#additive genetic variance is ~1
varp=1/h2
vare=varp-1
ngen=100

print("Base Population")
# a~ N(0,1)
a=randn(nqtl);

# p~ unif
p=rand(nqtl);
sum2pq=sum(2 .* p .* (1 .-p))
a=a/sqrt(sum2pq);

# storage for QTL genotypes
geno=Array{Int64,3}(undef,nanim,nqtl,2);
Z=Array{Int64,2}(undef,nanim,nqtl);

# simulate base animals
for i in 1:nanim
    for j in 1:nqtl
        for k in 1:2
            geno[i,j,k]=rand(Binomial(1,p[j]),1)[1]
        end
    end
end

vargen=Array{Float64,1}(undef,ngen)

for igen in 1:ngen

    println("-----------------------")
    println("generation number=", igen)
    println("-----------------------")
    for i in 1:nanim
        Z[i,:]=geno[i,:,1]+geno[i,:,2]
    end
    tbv=Z*a;
    println("var tbv= ",var(tbv))
    vargen[igen]=var(tbv)

    # generate phenotype with h2
    y=tbv+randn(nanim)*sqrt(vare);
    println("realized var y= ",var(y))
    println("realized h2= ",var(tbv)/var(y))

    # returns position in ranking based on y
    ranking=sortperm(y,rev=true);
    # pick best animals
    best=ranking[1:nselect];

    tbvselect=tbv[best];
    println("var tbvselect= ",var(tbvselect))

    # now create next generation sampling 2 alleles at random from 
    # selected individuals (possibly with selfing)
    genonew=Array{Int64,3}(undef,nanim,nqtl,2);
    for i in 1:nanim;
        for j in 1:nqtl
            for k in 1:2
                ancestor=rand(1:nselect);
                orig=rand(1:2);
                genonew[i,j,k]=geno[ancestor,j,orig];
            end
        end
    end
    geno[:,:,:]=genonew #not the same as geno=genonew !!!
end
plot(vargen)