Adrian became the third person to read my short story and - how shall we say? - not necessarily get the author's intent. The format - an epistolary structure based on blog entries - seems to me to be at fault and I think the answer is to remake the underlying story in a more conventional form.
Reading my forthcoming SMXR358 Summer School material, I came across the following method of estimating uncertainty when counting randomly occurring events (I quote).
If the number of randomly occurring events counted in a given period is N, the uncertainty in this count (i.e. the likely difference between N and the mean that would be found from a long series of repeated experiments) is √(N).
Example: if the mean is 100 counts, 68% 0f experiments will produce results within 90 and 110 counts i.e. +/- √100.
So where does that come from? I asked a tutor last year and we had an inconclusive discussion about Poisson distributions which seemed beside the point. Anyway, I reviewed it again and came to the following conclusion.
Suppose that there is an underlying binomial process with very low probability p of producing an observable event, and high probability q = (1-p)≈ 1 of not producing an observable event in a given time period. Radioactive decay would be an example.
Let T be the total number of underlying events (observable and not) which occur during the course of a single experiment and let Tp be the expected number of observable events. The average of the counts N which are measured would tend to this value over a number of experiments.
The standard deviation σ of the distribution is √(Tpq) = √(N *q) ≈ √N, as q is nearly 1.
So using the normal distribution approximation, 68% of measured values of N will lie within +/- 1σ of the mean, Tp, which is +/- √N as stated. Note that we use measured N as a best estimator for Tp.
I could have got to this more elegantly. The p ≈ 0 and q ≈ 1 condition and the large number of trials T is the condition for the Poisson distribution to apply. And in the Poisson distribution, the variance = the mean N, so the standard deviation is the square root of the mean = √N.