I'm making a game that has to simulate a lot of random events, but we only care about the sum total effect of those events. Calculating each event separately is costly. We can do better.

The naive approach

Imagine a sheep farm. On each tick in our simulation, there's a chance the sheep reproduce. We could do it like this:

local sheep_born = 0
for _ = 1, sheep_count do
  if math.random() < 0.05 then
    sheep_born = sheep_born + 1
  end
end

Thss is simple and it works, but it obviously starts to get very expensive as the number of sheep increases. Even more so if we have multiple farms.

The better approach

What if there was a magical function that could calculate the result of the above loop without calculating a random value thousands of times? I know practically nothing about statistics, so I was quite lost here. I was sure a solution exists, but what is it? Luckily, a friend pointed me towards the normal approximation of binomial distributions.

Implemented below in lua.

local function use_binomial_approximation(n, p)
  if n <= 5 then
    return false
  end
  if n >= 20 then
    -- Too many elements and the cost of precision isn't worth it
    return true
  end
  -- We allow a slightly larger skewness than usual
  local skew = math.abs(1 - 2 * p) / math.sqrt(n * p * (1 - p))
  return skew < 0.5
end

-- Produces a random number in the binomial sum distribution using
-- https://en.wikipedia.org/wiki/Binomial_distribution#Normal_approximation
-- Simulates:
--   x = 0
--   for i=1,n do
--     if math.random() < p then
--       x = x + 1
--     end
--   end
---@param n number number of uniform numbers to sum
---@param p number|nil chance each number is 1
---@return number
function math.binomial_distribution(n, p)
  p = p or 0.5

  if use_binomial_approximation(n, p) then
    local x = math.normal(n * p, n * p * (1 - p))

    -- This messes with the distribution slightly, but we want the result to
    -- never be outside the usual bounds
    x = math.clamp(x, 0, n)

    -- Round to nearest integer
    return math.floor(x + 0.5)
  else
    local x = 0
    for _ = 1, n do
      if math.random() < p then
        x = x + 1
      end
    end

    return x
  end
end

With this, we can convert the above snippet to

local sheep_born = math.binomial_distribution(sheep_count, 0.05)