Here’s a demo of the Bates distro.

library(tidyverse)
library(tictoc)

Iteratively take samples from the Uniform(0,1) distribution, drawing n_draws values each time and repeating this process for n_iterations. At each step, compute the mean across the draws obtained so far: begin with the uniform sample, then the mean of two uniform samples, then the mean of three, and so on – each time averaging across corresponding draws (1 to n_draws).

sim_bates_distro <- function(n_draws, n_iterations) {
  tibble(pos = 1:n_draws) |>
    mutate(draw = map(pos, ~ runif(n = n_iterations, min = 0, max = 1))) |>
    unnest(draw) |>
    group_by(pos) |>
    mutate(iteration = row_number(),
           cummean = accumulate(draw, ~ .x + .y) / iteration) |>
    ungroup() |>
    mutate(iteration_f = factor(iteration, labels = paste0("Step ", 1:n_iterations)))
}

Generate the data.

tic()
df_means <- sim_bates_distro(1e6, 6)
toc()
## 97.23 sec elapsed

This shows the first draw and cumulative averaging:

head(df_means, 6) |>
  select(pos, iteration, draw, cummean)

Make a picture.

ggplot(df_means, aes(x = cummean, fill = iteration)) +
  geom_histogram(bins = 500) +
  facet_wrap(~ iteration_f, scales = "free_y", ncol = 3) +
  labs(
    title = "An illustation of the the (Grace) Bates distribution",
    x = "",
    y = "Frequency"
  ) +
  theme(legend.position = "none")