Probability distribution describing Simpson's diversity

Question

Say we have $n$ species of unknown abundance in a community. Each species has the proportion $p_i$, such that $\sum_ip_i = 1$. Simpson's diversity is defined as $S = \sum_ip_i^2$, with a maximum value of 1 and a minimum value of $1/n$.

If we draw a sample from this population, I believe the proportions of each species we observe would be understood to be drawn from an $n$-dimensional Dirichlet distribution $D(p_1, p_2, ... p_n)$. Given this distribution, is there a way to derive a probability distribution for the diversity index $S$?

Answer 1

The problem of obtaining the distribution of the product (or sum) of Dirichlet distributions has been addressed multiple times(e.g., Nadarajah et. al. 2004). You can find more examples by looking for the distribution of the product and sum of Beta distributions, which is a particular case of the Dirichlet distribution. The distribution that results from this product and sum is commonly expressed in terms of Gaussian hypergeometric functions.

However, I found that Monte Carlo simulation provides a much simpler way of obtaining these distributions numerically. Consider the following example in R. In the particular case of two subpopulations, we can use the Beta distribution instead of the Dirichlet distribution:

f <- 0.2
n <- 20

reps <- 100000

S <- c()
for(r in 1:reps){
    p1 <- rbeta(1,f*n,(1-f)*n)  # Draw random values
    p2 <- rbeta(1,(1-f)*n,f*n)
    S[r] <- p1**2 + p2**2
}

hist(S,breaks=100)