How would one find the uncertainty in a mean if the data points themselves have zero-order uncertainty? Sorry if this question is this community's equivalent of asking a chef how to boil water, but if you had a data set that consists of:
[A±a, B±b, C±c, ..., N±n], where each value has a corresponding zero-order measurement uncertainty (as in, if you have a ruler with 1 mm ticks, the uncertainty would be 0.5 mm), you could very easily find the mean by just summing each nominal value and dividing by the number of values. But how would you find the uncertainty in that mean?
More specifically, say you had the following data set, obtained by repeated trials of a system:
[16±0.5, 21±0.5, 22±0.5]. The mean is obviously (16+21+22)/3 = 19.7. But how would you find the uncertainty of that mean?

Assuming the observations are collected independently of each other, the easiest way I can think of is to propagate uncertainty by using simulation. The idea is to generate random vectors from the (hyper-)cube and take the average of their coordinates; do this a large number of times and collect all the values obtained (the histogram below).
Here is an R code to illustrate it.
# set up a function to handle the simulation process
gen_x <- function() {
x = c(runif(1, 16-0.5,16+0.5),
      runif(1, 21-0.5,21+0.5),
      runif(1, 22-0.5, 22+0.5))
return(mean(x))
}

N=1e+4

# evaluate the function N times
hh <- sapply(1:N, function(x) gen_x())
hist(hh)
abline(v = mean(c(16,21,22)))


The histogram shows the distribution of the average of 16,21,22 taking into account the uncertainty of $\pm 0.5$ and assuming independent uniform distributions centred at the given values. The vertical line shows the sample average $(16+21+22)/3$.
The above solution gives the entire distribution of the sample average. If you are only interested in the variance of the sample average, then that's just equal to $1/(12n)$. Indeed, If you denote the sample $X_1,\ldots,X_n$, setting $\bar X = n^{-1}\sum_i X_i$ the sample average and assuming $X_i$ are independent, we have
$$
\text{var}(\bar X) = n^{-2}(\text{var}(X_1)+\cdots +\text{var}(X_n)) = n^{-2}\frac{n}{12} = (12n)^{-1}.
$$
Here I have used the fact that if $X \sim \text{Unif}(a,b)$, then $\text{var}(X) = (b-a)^2/12$.
So to sum up, if your $n$ samples are independent and the measures can be thought of as $\pm 0.5$ around the value measured, then, under the uniform assumption of the measure within each interval, the variance of the average is $(12n)^{-1}$.
