Phys 238:  Uncertainties resulting from Random, Independent errors

Clear["Global`*"]

Individual Random Numbers

RandomReal[{5 , 7}] (* Generates a random number between 5 and 7. *)

RandomReal[{5, 7}, 10] (* Generate a bunch (10) of them *)

vals = RandomReal[{5, 7}, 5000];  (* A big bunch *)

Histogram[vals]

Repeating a single measurement that results from the sum of many individual random numbers should give a normal distribution.

Define a function that represents one data point that is a result of many independent random numbers.

onepoint := Mean[RandomReal[{5, 7}, 200]]

This returns a different measurement each time you try it.

onepoint

onepoint

Now imaging doing 5000 such measurements.

data = Table[onepoint, {i, 1, 5000}];

Histogram[data]

μ = Mean[data]
σ = StandardDeviation[data]

How many are within ± 1 σ of the mean?  Use the ‘Select’ function to select data points within σ of μ.

Length[Select[data, μ - σ ≤ #  ≤ μ + σ &]]/Length[data];
StringForm["``, or ``%", %, 100 *N[%]]

Now try to model this distribution by a Gaussian.  Getting a manual list of the histogram numbers is a bit tedious.  The HistogramList[ ] function actually gives two separate lists, first is a list of bin boundaries, and second is the list of counts in each bin.

{bins, counts} =N[HistogramList[data]]

bins

counts

Interleaving them is straightforward, if verbose.  We will assign each height to the center point of each bin by constructing a new table.

hist = Table[{(bins[[i]] + bins[[i+1]])/2., counts[[i]]}, 
{i, 1, Length[counts]}];

Enter an unnormalized gaussian (Eq. 5.20 in Taylor’s text):

gauss[A_, μ_, σ_, x_] := (A/(σ Sqrt[2 π])) Exp[-((((x - μ))^2)/(2 (σ^2)))]

Clear[A, μ, σ]
(* The initial fit attempt does always work well.  It will be useful to seed the process with guesses at the mean and standard deviation. *)result = NonlinearModelFit[hist, gauss[A, μ, σ, x], 
{A, 
{μ, Mean[data]}, 
{σ, StandardDeviation[data]}
}, x]

result["BestFitParameters"]

result["ParameterConfidenceIntervalTable"]

These are very close to the calculated mean and standard deviation:

{Mean[data], StandardDeviation[data]}

Show the histogram, the data used in the fit, and the best-fit gaussian all on the same graph.  Add useful titles.

randerrplot = Show[{
Histogram[data],
ListPlot[hist, PlotStyle -> Black], Plot[result[x], {x, Min[bins], Max[bins]}]},
Frame->True, FrameLabel -> {"Value", "Count"}, LabelStyle -> Larger, ImageSize -> Scaled[0.8]]

SetDirectory[NotebookDirectory[]] (* Save it in this directory *)

Export["randerrplot.pdf", randerrplot]