Sample code 7.1. To generate repeatable random numbers, use SeedRandom[X], where X any integer.

For[i=1,i≤5,i++,
Print[i, " ",RandomReal[]]];(*end for i*)

For[i=1,i≤5,i++,
SeedRandom[2];
Print[i, " ",RandomReal[]]];(*end for i*)

Sample code 7.2. Random number is generated using RandomReal[] distributed uniformly in [0,1]. To generated 25 random number distributed according to Gaussian distribution NormalDistribution[μ,σ]
can be generated using RandomVariate[NormalDistribution[μ,σ],25].

data=RandomVariate[NormalDistribution[4,10],10];
Print[data]

Sample code 7.3. We can bin the 10^4 random number generated with RandomVariate[NormalDistribution[1,3],10^4] into 20 bins. This can be acchieved using the command Histogram[data,20,”ProbabilityDenstiy”].

data=RandomVariate[NormalDistribution[4,10],10^4];Histogram[data,20,"ProbabilityDensity"]

Sample code 7.3.1. Plot[PDF[NormalDistribution[4,10],x],{x,-40,40},PlotStyle->Thick] will plot the actual Gaussian Distribution with μ=4,σ=10 for x in [-40,40]
Overlapping the histogram and the Gaussian distribution curve, we can then check that indeed the random numbers are distributed according to the guassian distribution with the corresponding parameters  μ=4,σ=10.

Show[
Histogram[data,20,"ProbabilityDensity"],
Plot[PDF[NormalDistribution[4,10],x],{x,-40,40},PlotStyle->Thick]
](*end show*)

Sample code 7.4. 1-D random walk simulation involving only A SINGLE walker, simulated for a time span of nlast timestep. We denote the displacement of the walker at timestep n with repect to the initial displacement x0 (at x at n = 0) as x_n.

(*Begin of initialisation*)
x0=0;nlast=400;(*nlast = total time *)
x=x0;
(*End of initialisation*)

(*Setting the initial conditions for the walker at n=0*)
tnx={{0,x0}};tnxsq={{0,x0^2}};
(*Print[j," ",pairx[j]," ",tpairx[j]];*)
(*End of setting the initial conditions for the walker at n=0*)

For[n=1,n≤nlast,n++,
q=Random[Integer];
If[q==0,x=x+1,x=x-1];(*end if*)(* Shift the walker randomly*)
(*Print[q," ",x];*)
nx={n,x};nxsq={n,x^2};
tnx=Append[tnx,nx];tnxsq=Append[tnxsq,nxsq];
];(*end of n*)

ListPlot[{tnx(*,tnxsq*)},
(*Mesh->Automatic,*) PlotStyle->(*{PointSize[0.001],PointSize[0.003],PointSize[0.005],PointSize[0.007]},*)
{PointSize[0.005],Red,Blue,Yellow,Purple,Green},
PlotRange->All,PlotLabel->"Random walk of a single walker as a function of time"]
averx=Sum[tnx[[n]][[2]],{n,1,nlast}]/nlast;

Sample code 7.5. 1-D random walk simulation involving jlast walkers, simulated for a time span of nlast timestep.  Basically, code 7.5 and 7.4 are similar except that here the code evaluates jlast walkers instead of only one, as in code 7.4.  The code is obtained by modifying sample code 7.4 by adding an additional loop, indexed by j. We wish to know the following quantities at each time step n:
(i) xave_n = Sum[x_j,{j,1,jlast}]/jlast, at time n.
(ii) Average displacement squared at timestep n, xsqave_n=Sum[xsq_j,{j,1,jlast}]/jlast; 
Here, xave_n  is the average of all the displacements of the jlast walkers at time step n.

(*Begin of initialisation*)
x0=0;nlast=100;(*nlast = total time *)
jlast=500;(*jlast = total number of walker*)
(*End of initialisation*)

(*Setting the initial conditions for the walker at n=0*)
For[j=1,j≤jlast,j++,
x[j]=x0;
tnx[j]={{0,x0}};tnxsq[j]={{0,x0^2}};
];(*end for j*)
tnxave={{0,0}};tnxsqave={{0,0}};tnvarx={{0,0}};
(*End of setting the initial conditions for the walker at n=0*)


For[n=1,n≤nlast,n++,
For[j=1,j≤jlast,j++,
q=Random[Integer];
If[q==0,x[j]=x[j]+1];
If[q==1,x[j]=x[j]-1];(*end if*)(* Shift the walker randomly*)
nx[j]={n,x[j]};nxsq[j]={n,x[j]*x[j]};
tnx[j]=Append[tnx[j],nx[j]];tnxsq[j]=Append[tnxsq[j],nxsq[j]];
];(*end for j*)
avex=Sum[x[j],{j,1,jlast}]/jlast;
avexsq=Sum[x[j]^2,{j,1,jlast}]/jlast;
tnxave=Append[tnxave,{n,avex}]; 
tnxsqave=Append[tnxsqave,{n,avexsq}]; 
varx=avexsq-avex^2;
tnvarx=Append[tnvarx,{n,varx}]; 
];(*end of n*)

After the random walk of jlast walkers been simulated, their averaged displacement, average displacement squared and variance as functions of time can be displayed using ListPlot as follows.

j1=1;j2=jlast;
ListPlot[
Table[tnxsq[j],{j,j1,j2}],
(*Mesh->Automatic,*) PlotStyle->PointSize[0.01],
PlotRange->All,PlotLabel->"xsq as a function of time for all jlast walkers. At a given n, the averaged value of xsq is just given by Sum[xsq[j],{j,1,jlast}]"];

ListPlot[
{tnxave},Joined->True,
Mesh->Automatic, PlotStyle->{{PointSize[0.01],Red}},
PlotRange->All,PlotLabel->"{x averaged over all walkers as a function of time}"];

ListPlot[
{tnxsqave},Joined->True,
Mesh->Automatic, PlotStyle->{{PointSize[0.01],Red}},
PlotRange->All,PlotLabel->"xsq averaged over all walkers as a function of time"];

ListPlot[
{tnvarx},Joined->True,
Mesh->Automatic, PlotStyle->{{PointSize[0.01],Red}},
PlotRange->All,PlotLabel->"variance of x as a function of time"]

A detail analysis of the random walk on the jlast 1-D walkers shows that 
<x^2>_n = 2Dn, where D is a constant called the diffusion coeffience, n denote time (in unit of multiple of 1). For 1D case, D = 1/2, so that <x^2>_n = n, which can be compared easily with the simulation result.

Dcoef=1/2;
y[n_]:=2Dcoef* n;Plot[y[n],{n,0,nlast}];ytheory=Table[{n,y[n]},{n,0,nlast}];

ListPlot[
{ytheory,tnvarx(*tnxsqave*)},
(*Mesh->Automatic,*) PlotStyle->{{PointSize[0.01],Red},{PointSize[0.01],Green}},
PlotRange->All]

In the above code, we also have simulated the variance of the walker, defined as . We can plot variance as a function of n. According to analysis on the 1-D jlast random walker, viriance will depend on time via the relation =, where l = 1, the spatial step the walkers are shifted in each simulation step.

y[n_]:= n;Plot[y[n],{n,0,nlast}];ytheory=Table[{n,y[n]},{n,0,nlast}];
ListPlot[
{ytheory,tnvarx},
(*Mesh->Automatic,*) PlotStyle->{{PointSize[0.01],Red},{PointSize[0.01],Green}},
PlotRange->All]

Sample code 7.6. 2-D random walk simulation involving jlast walkers, simulated for a time span of nlast timestep. We would like to visual the walkers on a 2D plan, using the Manipulate function of mathematica. We assume all walkers are located at r0={x0,y0} at n=0.

(*Begin of initialisation*)
x0=0;y0=0;
nlast=100;(*nlast = total time *)
jlast=500;(*jlast = total number of walker*)
(*End of initialisation*)

(*Setting the initial conditions for the walker at n=0*)
For[j=1,j≤jlast,j++,
x[j]=x0;y[j]=y0;r[j]={x0,y0};
tnx[j]={{0,x0}};tnxsq[j]={{0,x0^2}};
tny[j]={{0,y0}};tnysq[j]={{0,y0^2}};
trn[j]={{0,0}};trnsq[j]={{0,0}};
tn[j]={{0,0}};
];(*end for j*)

For[n=1,j≤nlast,n++,
tr[n]={{0,0}};
];(*end for n*)

tnxave={{0,0}};tnxsqave={{0,0}};tnvarx={{0,0}};
tnyave={{0,0}};tnysqave={{0,0}};tnvary={{0,0}};
tnrave={{0,0}};tnrsqave={{0,0}};tnvarr={{0,0}};

(*End of setting the initial conditions for the walker at n=0*)


For[n=1,n≤nlast,n++,
For[j=1,j≤jlast,j++,
q=Random[];
If[q≤0.25,x[j]=x[j]+1;y[j]=y[j]+1];
If[0.25<q≤0.5,x[j]=x[j]-1;y[j]=y[j]+1];
If[0.5<q≤0.75,x[j]=x[j]+1;y[j]=y[j]-1];
If[0.75<q≤1.00,x[j]=x[j]-1;y[j]=y[j]-1];
nx[j]={n,x[j]};nxsq[j]={n,x[j]^2};
ny[j]={n,y[j]};nysq[j]={n,y[j]^2};

r[j]=Sqrt[x[j]^2+y[j]^2];
tnr[j]=Append[tnr[j],r[j]];tnrsq[j]=Append[tnrsq[j],r[j]];
];(*end for j*)
avex=Sum[x[j],{j,1,jlast}]/jlast;
avexsq=Sum[x[j]^2,{j,1,jlast}]/jlast;
tnxave=Append[tnxave,{n,avex}]; 
tnxsqave=Append[tnxsqave,{n,avexsq}]; 
varx=avexsq-avex^2;
tnvarx=Append[tnvarx,{n,varx}]; 

avey=Sum[y[j],{j,1,jlast}]/jlast;
aveysq=Sum[y[j]^2,{j,1,jlast}]/jlast;
tnyave=Append[tnyave,{n,avey}]; 
tnysqave=Append[tnysqave,{n,aveysq}]; 
vary=aveysq-avey^2;
tnvary=Append[tnvary,{n,vary}]; 

aver=Sum[r[j],{j,1,jlast}]/jlast;
aversq=Sum[r[j]^2,{j,1,jlast}]/jlast;
tnrave=Append[tnrave,{n,aver}]; 
tnrsqave=Append[tnrsqave,{n,aversq}]; 
varr=aversq-aver^2;
tnvarr=Append[tnvarr,{n,varr}]; 

tr[n]=Table[{x[j],y[j]},{j,1,jlast}];
];(*end of n*)

After the 2D random walk is done, we have collected all the various value of r[j] for all n. Based on these values we can caculate easily <r^2>, which is stored in tnrsqave.

Manipulate[
ListPlot[
tr[ndum],Frame->True,PlotStyle->{PointSize->0.005},PlotRange->{{-75,75},{-75,75}}](*end graphics*)
,{ndum,1,nlast,1}](*end manipulate*)

After the 2D random walk is done, we have collected all the various value of r[j] for all n. Based on these values we can caculate easily <r^2>, which is stored in tnrsqave. The graph of <r^2> is plotted against n to show a linear relation with n. It is easily found that the slop for the 2D random walk in <r^2> vs. n is given by <r^2> = 2*D*n, where D = 1 in this case.

y[Dr_,n_]:=2*Dr*n;
fit[Ddum_]:=Table[{n,y[Ddum,n]},{n,1,nlast}];
ListPlot[
{fit[1],tnrsqave},PlotLabel->"<r^2> vs. n. Dr = 1"
]

The variance in r is also plotted using the data stored in tnvarr. We find that D = 0.20.

y[Dr_,n_]:=2*Dr*n;
fit[Ddum_]:=Table[{n,y[Ddum,n]},{n,1,nlast}];
ListPlot[
{fit[0.20],tnvarr},PlotLabel->"varr vs. n. Dr = 0.20"
]

y[Dx_,n_]:=2*Dx*n;
fit[Ddum_]:=Table[{n,y[Ddum,n]},{n,1,nlast}];
ListPlot[
{fit[1/2],tnxsqave},PlotLabel->"<x^2> vs. n. Dx = 1/2"
]

y[Dy_,n_]:=2*Dy*n;
fit[Ddum_]:=Table[{n,y[Ddum,n]},{n,1,nlast}];
ListPlot[
{fit[1/2],tnysqave},PlotLabel->"<y^2> vs. n. Dy = 1/2"
]