How to sample from a distribution with discrete variables with a known average?

Question

My motivation is to produce carcasses of animals in an ecosystem.

• The animals have discrete sizes in kg (75, 216, 700, 2500, 5000, 8500, 25000).
• I also have the estimated percentage each animal contributes to the total ecosystem based on scaling relationships (49.3, 36.8,6,6.7,0.6,0.4,0.2)
• Finally, I also have the average amount of kg the system should support (1752kg)

So I want to have a discrete distribution that uses the first two bits of data to produce a system with 1752kg of animals on average. So it's fine that sometimes it would overshoot or undershoot the average.

This answer gets some of the way there but my average value seems to make a difference.

I ultimately want to implement this in NetLogo code but I'm familiar with R if people want to illustrate their answers that way.

Answer 1

You don't have much details about the process, so you can aim only at a rough approximation of the distribution. The simple, approximate solution is not hard though.

You know how often animals of different sizes appear in the population, we need those frequencies to hold in the samples. There is another constraint for the average weight of all the animals in the sample to be 1752 kg.

Let's denote the weight per animal type as $x_i$ and the "estimated percentage each animal contributes to the total ecosystem" to be $w_i$, then it's a categorical distribution with $\Pr(X = x_i) = w_i$. From here, we can find the expected weight of a single animal $E[X] = \sum_i x_i w_i = 439.96$. We can also calculate how many animals we need to have so that they on average have the desired weight. You want to sample $n$ animals such that their total weight is $T = 1752$, so $E[n X] = T$, and by the linearity of expectation, $E[nX] = nE[X]$, so we can find $T / E[X] = n$ and it is equal to $1752 / 439.96 = 3.98$. So you need to sample with replacement $\approx4$ animals with the probabilities $w_i$.

In R code, this is

> x <- c(75, 216, 700, 2500, 5000, 8500, 25000)
> w <- c(49.3, 36.8, 6, 6.7, 0.6, 0.4, 0.2)
> w <- w/sum(w)
> summary(replicate(100000, sum(sample(x, size=4, prob=w, replace=TRUE))))
##   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    300     582     723    1760    2725   50291 

Notice that the mean is $1760$ what is consistent with $4 \times E[X]$.

> sum(x * w)
## [1] 439.963
> 4 * sum(x * w)
## [1] 1759.852

This is not equal to $1752$ what shows us that either the numbers you have are not exact, or the system is more complicated and you would need to make more assumptions for a better approximation.

One thing you could do is to have $n$ random instead of fixed. For example, you could assume that $n$ follows a Poisson distribution with $\lambda = 3.98$, to get a much better approximation:

> Ex <- sum(x * w)
> T <- 1752
> En <- T/Ex
> summary(replicate(100000, sum(sample(x, size=rpois(1, En), prob=w, replace=TRUE))))
##   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      0     366     775    1752    2650   51907

Notice that we got better approximation for the mean, but at the cost of making an assumption about the distribution of the counts. This assumption may or may not have sense for the problem. The same applies for any other approach to the simulation that you would take: you have little details, so you would need to make smaller or larger assumptions about the process and do it wisely.