Randomness

AP1601
M. E. Mauel

Introduction

In our study of magnetism and protein folding, we will be using Random[] within our Monte Carlo simulation of thermaldynamic fluctuations in a solid.

In this notebook, we briefly evaluate the use of random numbers. We will illustrate the well-known rule that the uncertainty of samples decreases as the square-root of the number of samples.

For example, in our Ising simulations with a N × N lattice, the uncertainty of any average scales as 1/N. Averages of a 10 × 10 lattice will be accurate to 10%; averages of a 100 × 100 grid will be accurate to ± 1%.

Moments of Random Numbers

One of the most used tests of random numbers is to take moments. For uniformly distributed random numbers between 0 and 1, the kth moment is defined as 
	(1/N)(∑_(i=1)^N)(x^k) 
which equals 1/(k + 1) when N -> ∞.

nTable[n_Integer]:=Table[RandomReal[],{i,1,n},{j,1,n}]

t1=nTable[200];
ListDensityPlot[t1,Mesh->False,ColorFunction->Hue]

moment[t_List,k_Integer] := (Plus @@ Flatten[t^k])/Length[Flatten[t]]

p1=ListPlot[Table[{k,moment[t1,k]},{k,1,25}],PlotRange->{0.,0.5},PlotStyle->{RGBColor[0,0,1],PointSize[0.02]}];
p2=Plot[(1/(k+1)),{k,1,25}];Show[p1,p2,PlotLabel->"Average Moment of Order k"]

Changing the size of the grid…

p1=ListPlot[Table[{n,Abs[moment[nTable[n],3]-(1/(3+1))]},{n,20,200}],PlotRange->All,PlotStyle->{RGBColor[0,0,1],PointSize[0.02]}];
p2=Plot[(1/n),{n,20,200}];Show[p1,p2,PlotLabel->"Variation of Average Moment vs. N"]

Standard Deviation (or Variance)

The standard measure of fluctuations is the standard deviation or variance. The variance, (((Δx))^2), is defined as the square of the standard deviation, Δx.
	 (((Δx))^2) = ⟨(x^2)⟩ - ((⟨x⟩)^2)
where ⟨…⟩ means an average over the samples.

For an N × N lattice of random numbers, the standard deviation will be nearly constant.  As N -> ∞, Δx -> (Sqrt[3])/ 6 ~ 0.288675.

We can write a simple expression to calculate the variance, or we can use the standard package "StatisticsDescriptiveStatistics".

stddev[t_List] := Sqrt[moment[t,2] - moment[t,1]^2]

ListPlot[Table[{n, stddev[nTable[n]]}, {n, 20, 200}], PlotRange->{0.28,0.30}, PlotStyle->{RGBColor[0,0,1],PointSize[0.02]},PlotLabel-> "Std Deviation vs. N"]

The built-in functions in the Mathematica execute more quickly…

ListPlot[Table[{n, StandardDeviation[Flatten[nTable[n]]]}, {n, 20, 200}], PlotRange->{0.28,0.30}, PlotStyle->{RGBColor[0,0,1],PointSize[0.02]},PlotLabel-> "Std Deviation vs. N"]

Correlations

For true random numbers, there is no average correlation between any (different) pair of random numbers with zero mean.

We define a correlation between two random numbers as
	C(k) = (1/N)(∑_(i=1)^N)(Subscript[x, i])(Subscript[x, i+k])
If the numbers, (Subscript[x, i]), are truly random and uniformly distributed between 0 and 1, then C(k) = (1/4)for all k.

correlate[t_List,k_Integer] := (Plus @@ Flatten[t RotateRight[t,{k,0}]])/Length[Flatten[t]]

ListPlot[Table[{k,correlate[t1,k]},{k,1,50}],PlotRange->{0.,0.5},PlotStyle->{RGBColor[0,0,1],PointSize[0.02]}]

ListPlot[Table[{k,(correlate[t1,k] - 0.25) 100}, {k, 1, 50}],PlotStyle->{Red,PointSize[0.02]}, PlotLabel-> "Difference between ideal (%)"]

Very long lists of Random Numbers

In order to look at very long lists of random numbers, I choose not to store the numbers in a long list. Instead, we compute the standard deviation in a long loop. This requires processing power, but little memory.

myStdDev[n_Integer]:=Module[{x,sum1=0.,sum2=0.},Do[x=RandomReal[];sum1+=x;sum2+=(x^2),{n}];(Sqrt[(sum2/n)-((((sum1/n)))^2)])]

myStdDev[10^6] // Timing

Using Compile[…], we can significantly increase computation speed…

fastStdDev=Compile[{{n,_Integer}},Module[{x,sum1=0.,sum2=0.},Do[x=RandomReal[];sum1+=x;sum2+=(x^2),{n}];(Sqrt[(sum2/n)-((((sum1/n)))^2)])]]

fastStdDev[10^6] // Timing

diffData = Table[{(10^n),Abs[fastStdDev[(10^n)]-(Sqrt[(1/3)-(1/4)])]},{n,1,7, Log[10,2]}];
ListLogLogPlot[diffData,Joined-> True];
ListLogLogPlot[diffData,PlotStyle->{Red,PointSize[0.02]}];
dataPlt = Show[%%,%,PlotLabel->"Std Deviation Diff vs. N", AxesLabel -> {"N","Diff"}]

Fit[diffData, {1, 1/Sqrt[n]}, n]
sqCoeff = Coefficient[%, 1/Sqrt[n]]

fitData[n_] := sqCoeff/Sqrt[n]

LogLogPlot[fitData[n], {n, 10, (10^7)}];
Show[dataPlt, %]

Summary

We will use two common statistical tools in our statistical analysis of the Ising lattice.

First, we will take averages. The accuracy of this process is of order 1/N where N is the dimension of one side of the square lattice. (Accuracy scales like 1/(Sqrt[N×N]), the inverse square root of the total number of cells.)

Second, we will find two-point correlations. For random numbers, this correlation is independent of the separation between points and has a value of 0.25. Since a correlation is like an average, accuracy also requires a large number of lattice sites.