NeiH(group1)
NeiH(group2)

B <- 5000 # number of bootstrap replicates

# Initialization of three vectors that will receive the statistics for each bootstrap sample:
g1_boot <- rep(NA, B); g2_boot <- rep(NA, B); diff_boot <- rep(NA, B)

# Resampling (with replacement) in each population:
for (i in 1:B) {
    g1star <- sample(group1, replace=TRUE, size=length(group1)) # draw a random sample in group1
    g2star <- sample(group2, replace=TRUE, size=length(group2))
    g1_boot[i] <- NeiH(g1star)
    g2_boot[i] <- NeiH(g2star)
    diff_boot[i] <- g1_boot[i] - g2_boot[i]
}

# Histogram for the difference between the two indicators of gene diversity:
hist(diff_boot, xlim=c(-0.5,0.1), main="Difference between Nei's H indices")
# Bootstrap 95% IC for the difference:
quantile(diff_boot, 0.025)
quantile(diff_boot, 0.975)

hist(g2_boot)
abline(v=0.97, lty=2, lwd=2)