generateMM[] := SetAccuracy[RandomReal[NormalDistribution[0.53,0.014]],4]

How to describe and analyze distributions :
  A Mathematica walkthrough

Let's set up an engineering problem.    You are working in the M&M candy factory.  Your boss has asked you to set the machine to manufacture exactly 10,000 M&Ms, and he wants to know how much chocolate the machine will consume.   He doesn't want to over-buy chocolate (the excess would go to waste), nor underbuy (the company would have to place an expensive rush-order to finish the production).    "In fact", says the boss, "tell me the smallest amount of chocolate to buy such that I have only a 16% chance of running out.   Here are 100 M&Ms to measure." 
  
So, you need your best guess at the *average* amount of chocolate in 10, 000 M & Ms, and the error on that guess.    Let's start by measuring our 100 M&Ms.   I have a mathematica function that will (randomly) tell you the mass of one M&M, in grams, as measured on a scale with    Here it is:

generateMM[]

dataset = Table[generateMM[],{i,1,100}]

Histogram[dataset]

OK, so it looks like the average M & M is something like 0.52 or 0.54 grams.   But they aren't all exactly 0.53 or 0.54-- - there's a range.   This is important.    Let's check the average, range, etc.:

xmean = Mean[dataset]
xsigma = StandardDeviation[dataset]
xrange = {Max[dataset], Min[dataset]}

The first thing to notice is : the range is much greater than the precision of the scale.    If it' s really a 1 - mg - precision scale, then this is telling us that the M & Ms actually vary individually in mass.    Or, it could be telling us that the scale has a ~13 milligram error bar we didn' t know about.   (Can you think of a way to tell the difference?)

Anyway : Let' s tell the boss how much chocolate to buy to make 100 M & Ms.   Maybe we just say : these M & Ms added up to (x) and tell him to buy 100 times that.     But this doesn' t tell him the odds of over or underbuying.

guess1 = 100*Total[dataset]

Maybe we say : buy 10, 000 x the average of this set.   (Oh, wait, that' s the same thing, by definition).

guess2 = 10000* Mean[dataset]

Here' s the problem.That is the mean of this particular set of 100 M & Ms.These are probably more or less typical of any 100 M & Ms, but they' re not exactly precisely the same.The next 100 M & Ms will have a different mean, and the next will be different again.Your 100 - point mean is an estimator of the "true mean", but there' s an error on your estimate.

A common mistake is to look at the distribution of single M&Ms within your dataset, and suppose that this is the error on your knowledge of the mean.     Let's follow this (erroneous) line of thought for a moment and see where it goes wrong.  

Here's the bad approach in action:

Print["Mean = ", 10000*xmean, "+/-", 10000*xsigma , " grams"]
Print["The factory needs to buy ", 10000*(xmean+ xsigma), " grams"]
Print["and therefore has a 16% chance of running out." ]

Let' s  simulate a bunch of 10, 000 M & M run and see how often the factory runs out.   Here is one such run, which will output "0" if the factory has run out based on your estimate, "1" if it is OK (overbought)

diditfail[] := If[Total[ Table[generateMM[],{i,1,10000}]] > 10000*(xmean+xsigma),0,1]
diditfail[]

Let' s repeat the whole experiment (draw 100, measure mean, predict whole run, check for shortfall)  100 times and look for the purported 16 % chance of failure :

Table[
trialdataset = Table[generateMM[],{i,1,100}];
trialxmean = Mean[trialdataset];
trialxsigma = StandardDeviation[trialdataset];
diditfail[] := If[Total[ Table[generateMM[],{i,1,10000}]] > 10000*(trialxmean+trialxsigma),0,1];
diditfail[],
{i,1,100}]

Uh - oh.   You' re fired.  You overestimated the error, and your boss just bought (and threw away) excess chocolate with what looks like 100 % probability.     
   
   What went wrong?    Let's understand.   Look at your dataset again:

histo = Histogram[dataset]

xmean = Mean[dataset]
xsigma = StandardDeviation[dataset]
Print["I claim that the mean = ", xmean, "+/-", xsigma , " grams"]

That' s the same as the following statement :

Print["This is 100 datapoints drawn from a distribution whose true center is (with 68% probability) somewhere between ", xmean+xsigma, " and ", xmean-xsigma]

That' s a strange statement.   Let' s look at a "true" distribution whose mean is actually at the 68 % - likely - lower - bound, which the above statement says is a reasonably likely model for the data.

proposedmean = xmean-xsigma
Show[histo,Plot[25*E^(-((x-proposedmean)/xsigma)^2),{x,0.45,0.6}]]

If the true mean is really that low, it' s extraordinarily unlikely that so many data points would be that high.   A data/theory disagreement this bad does NOT happen 68 % of the time.   

What's going on?   THERE ARE TWO SEPARATE QUANTITIES that you must not conflate.  One is "sigma", the width of your distribution.   The other is the "error bar", the uncertainty in your knowledge of something.    

In this notebook, xsigma is the unpredictability of the mass of any individual M&M.     If you have measured only one M&M's mass, if is reasonably consistent with the idea that the mean is that-mass plus or minus sigma (assuming you knew sigma from somewhere else).   For example, in the plot below, the *single* datapoint is a reasonably plausible result of *any* of the three distributions suggested, centered as far as xsigma away from the mean:

verysmalldataset= {generateMM[]}
smallmean = Total[verysmalldataset];
Show[Histogram[verysmalldataset,{0.25,0.75,0.025}],Plot[2*E^(-((x-smallmean+xsigma)/xsigma)^2),{x,0.25,0.75},PlotRange->All,PlotStyle->Red],
Plot[2*E^(-((x-smallmean-xsigma)/xsigma)^2),{x,0.25,0.75},PlotRange->All,PlotStyle->Blue],
  Plot[2*E^(-((x-smallmean)/xsigma)^2),{x,0.25,0.75},PlotRange->All]
]

So, here is a fact I' ll state, but not derive.   YOUR KNOWLEDGE OF THE MEAN GETS BETTER AND BETTER  THE MORE DATAPOINTS YOU'VE USED TO DERIVE IT.   YOUR ERROR BARS ARE SUPPOSED TO BE STATEMENTS ABOUT THIS KNOWLEDGE, NOT JUST ABOUT THE WIDTHS OF DISTRIBUTIONS OF DATA POINTS.    The specific form of the "getting better" is as follows.   The error estimate is called the "standard error on the mean", and it's different than the "standard deviation" which is just the width of the distribution.

numberofdatapoints = 100
errorestimate = xsigma/Sqrt[numberofdatapoints]

Let' s repeat the ugly - looking histogram above and see how much better this looks.  This is a distribution which is plausibly (at 68 % level) the source of these 100 data points, even though the distribution is centered one - error - bar too low.

proposedmean = xmean-errorestimate
Show[histo,
Plot[35*E^(-((x-proposedmean)/xsigma)^2),{x,0.45,0.6}]]

So, let' s go to the boss of the M & M factory with our actual error estimate.

Print["The mean M&M weight is = ", xmean, "+/-", xsigma/Sqrt[100] , " grams"]
Print["Individual M&Ms vary in mass by +/- " ,xsigma]
Print["The factory needs to buy ", 10000*(xmean+xsigma/Sqrt[100]  ), " grams of chocolate"]
Print["in order to have a 16% chance of running out." ]

Table[
trialdataset = Table[generateMM[],{i,1,100}];
trialxmean = Mean[trialdataset];
trialxsigma = StandardDeviation[trialdataset];
diditfail[] := If[Total[ Table[generateMM[],{i,1,10000}]] > 10000*(trialxmean+trialxsigma/Sqrt[100]),0,1];
diditfail[],
{i,1,100}]

Yes, as desired, that was a conservative - enough chocolate purchase that the runs succeeded about 82 % of the time.  That is what error bars are supposed to mean.    

Some notes: 

a) Yes, you can derive this 1/sqrt(N) behavior from the physics.   When you're taking the mean of three data points---let's call them a,b,c, for which N=3---you are computing

m = 1/N * (a+b+c)

Suppose that a has an error of D, as does b, as does c.   Do error propagation:

dm^2 = ((dm/da)*D)^2 + ((dm/db)*D)^2 + ((dm/dc)*D)^2) 

and you get that dm = D/sqrt(N).    The error/spread on the mean is smaller than the error/spread on the individual datapoints you're using to calculate it.   


b) In this case, is the analysis different if (for example) the error were due to the scale, rather than due to the M&Ms varying individually?   Not really.   The fact is the *any scatter* makes it hard to estimate where the mean is, and that's true no matter where the scatter is coming from.

c) In this case, you don't know where the "true" mean of your 100 M&Ms is.   However, even if you knew the true mean, on any given 10,000-candy factory run, you'd get statistical fluctuations in the total of the 10,000.    Is that a problem?    Yes, you're very smart.    "manufacuring 10,000 M&Ms" is a statistical process just like "picking 100 to measure" is.    However, the error on the manufacturing will drop as 1/sqrt(10,000), which is a good but smaller than 1/100, so in this example it really is the measurer's inability to estimate the mean that dominates.    

d) What labs should pay attention to this notebook?   ALL OF THEM.   Anything where you're repeating a measurement, or combining multiple similar measurements, hoping that the "average" is a better answer than "the first run".    It probably is, but only if you correctly use the average to reduce your error bar.