Clear all symbols from previous evaluations to avoid problems

Clear["Global`*"]

Generating random variables with a specified distribution

Special methods for the normal distribution

The normal distribution N(μ,(σ^2))

ρ(x)=(1/(σ(Sqrt[2π])))exp(-(((x-μ))^2)/2(σ^2))

can be efficiently generated for example using the central limit theorem. If we sum up 12 uniformly distributed (μ=1/2, (σ^2)=1/12) random variables we get N(6,1). Thus by shifting to 0 and then by multiplying by σ and adding μ we get the general normal distribution:

myNormalSums[n_,par_]:=Module[{out=ConstantArray[0.0,{n}]},
Do[
out[[i]]=0.0;
Do[
out[[i]]=out[[i]]+RandomReal[],
{j,1,12}];
out[[i]]=par[[1]]+par[[2]](out[[i]]-6.0),
{i,1,n}
];
Return[out]
];

μ=10; σ=3;
n=10^5;
points=myNormalSums[n,{μ,σ}];
Show[Histogram[points,{0,20,1/4},"ProbabilityDensity"],Plot[PDF[NormalDistribution[μ,σ],x],{x,0,20}]]
expectedValue=Total[points]/n ;
variance=Total[points^2]/n-expectedValue^2 ;
Print["Expected value and standard deviation should be: ", N[μ],", ",N[σ]];
Print["Expected value and standard deviation obtained:  ", expectedValue,", ",Sqrt[variance]];
kmax=10; corr=ConstantArray[1.0,kmax+1];
Do[
corr[[k+1]]={k,(points[[1;;n-k]].points[[k+1;;n]]/(n-k)-expectedValue^2)/variance},
{k,0,kmax}
];
ListPlot[corr,PlotRange->All]

Another efficient way of generating the normal distribution is to consider a Gaussian distribution in two dimensions and use the polar coordinates to transform it to a product of exponential and uniform distribution:

exp(-((Subscript[x, 1]^2)+(Subscript[x, 2]^2))/2)d(Subscript[x, 1])d(Subscript[x, 2])∼exp(-(r^2)/2)rdrdθ∼exp(-u)dudθ

If u is generated between 0 and ∞ with an exponential distribution using

y(u)=(∫_0^u)exp(-t)dt=1-exp(-u)    ⟹    u(y)=-ln(1-y)

and θ is uniformly distributed between 0 and 2π then

(Subscript[x, 1])=(Sqrt[2u])cos(θ), (Subscript[x, 2])=(Sqrt[2u])sin(θ)

are distributed normally.

myNormal2D[n_,par_]:=Module[{r,θ,out=ConstantArray[0.0,{n+Mod[n,2]}]},
Do[
r=Sqrt[-2Log[1-RandomReal[]]];
θ=2π RandomReal[];
out[[i]]=par[[1]]+par[[2]]r Cos[θ];
out[[i+1]]=par[[1]]+par[[2]]r Sin[θ],
{i,1,n+Mod[n,2],2}
];
Return[out[[1;;n]]]
];

μ=10; σ=3;
n=10^5;
points=myNormalSums[n,{μ,σ}];
Show[Histogram[points,{0,20,1/4},"ProbabilityDensity"],Plot[PDF[NormalDistribution[μ,σ],x],{x,0,20}]]
expectedValue=Total[points]/n ;
variance=Total[points^2]/n-expectedValue^2 ;
Print["Expected value and standard deviation should be: ", N[μ],", ",N[σ]];
Print["Expected value and standard deviation obtained:  ", expectedValue,", ",Sqrt[variance]];
kmax=10; corr=ConstantArray[1.0,kmax+1];
Do[
corr[[k+1]]={k,(points[[1;;n-k]].points[[k+1;;n]]/(n-k)-expectedValue^2)/variance},
{k,0,kmax}
];
ListPlot[corr,PlotRange->All]

von Neumann method

myRandomNeumann[ρ_,a_,b_,R_,n_]:=Module[{g,x,y,rejected=0,out},
out=ConstantArray[0.0,n];
Do[
g=1;
While[g==1,
x=a+RandomReal[](b-a);
y=R RandomReal[];
If[y<ρ[x],g=0,rejected++]
];
out[[i]]=x,
{i,1,n}
];
Print[rejected," rejected points to generate ", n, " points."];
Return[out]
];

Normal distribution generated by von Neumann method:

a=0;b=20;
μ=10; σ2=3;
ρ:=PDF[NormalDistribution[μ,σ2],#]&;
R:=PDF[NormalDistribution[μ,σ2],μ];
Plot[{ρ[x],R},{x,a,b}]
n=10^4;
points=myRandomNeumann[ρ,a,b,R,n];
Show[Histogram[points,{a,b,1/4},"ProbabilityDensity"],Plot[PDF[NormalDistribution[μ,σ2],x],{x,a,b},PlotStyle->Red]]
expectedValue=Total[points]/n ;
variance=Total[points^2]/n-expectedValue^2 ;
Print["Expected value and standard deviation should be: ", N[μ],", ",N[σ]];
Print["Expected value and standard deviation obtained:  ", expectedValue,", ",Sqrt[variance]];
kmax=10; corr=ConstantArray[1.0,kmax+1];
Do[
corr[[k+1]]={k,(points[[1;;n-k]].points[[k+1;;n]]/(n-k)-expectedValue^2)/variance},
{k,0,kmax}
];
ListPlot[corr,PlotRange->All]

Exponential distribution generated by von Neumann method (not very efficient):

a=0;b=10;
ρ:=Exp[-#]&;
R:=1;
Plot[{ρ[x],R},{x,a,b}]
n=10^4;
points=myRandomNeumann[ρ,a,b,R,n];
Show[Histogram[points,{a,b,1/4},"ProbabilityDensity"],Plot[Exp[-x],{x,a,b},PlotRange->{0,1},PlotStyle->Red]]
expectedValue=Total[points]/n ;
variance=Total[points^2]/n-expectedValue^2 ;
Print["Expected value and standard deviation obtained:  ", expectedValue,", ",Sqrt[variance]];
kmax=10; corr=ConstantArray[1.0,kmax+1];
Do[
corr[[k+1]]={k,(points[[1;;n-k]].points[[k+1;;n]]/(n-k)-expectedValue^2)/variance},
{k,0,kmax}
];
ListPlot[corr,PlotRange->All]

Metropolis algorithm

One dimensional case, ρ must be a function of x:

myRandomMetropolis1D[n_,ρ_,a_,b_,x0_,δ_]:=
Module[{x=x0,xtrial,ratio,ρx,out=ConstantArray[0.0,n],accepted=0},
ρx=ρ[x];
Do[
(* try a new random step *)
xtrial=x+δ RandomReal[{-1,1}];
(* check if we are still in the given region, if not, adjust *)
If[xtrial<a,xtrial=a+δ RandomReal[]];
If[xtrial>b,xtrial=b-δ RandomReal[]];
ratio=ρ[xtrial]/ρx;
(* decide whether to make this step according to the distribution *)
If[ratio⩾1||ratio>RandomReal[],x=xtrial;ρx=ρx ratio;accepted++];
out[[i]]=x,
{i,1,n}
];
Print["Number of accepted trials: ",accepted, " from ", n];
Return[out]
];

Applied to the normal distribution:

n=10^5;
a=0;b=20;
μ=10; σ=3;
ρ:=PDF[NormalDistribution[μ,σ],#]&;
x0=10.0; step=5.0;
points=myRandomMetropolis1D[n,ρ,a,b,x0,step];
data=RandomReal[NormalDistribution[μ,σ],{n}];
ListPlot[{data,points},Joined->{False,True}]
Show[Histogram[points,{a,b,1/4},"ProbabilityDensity"],Plot[PDF[NormalDistribution[μ,σ],x],{x,a,b},PlotStyle->Red]]
expectedValue=Total[points]/n ;
variance=Total[points^2]/n-expectedValue^2 ;
Print["Expected value and standard deviation should be: ", N[μ],", ",N[σ]];
Print["Expected value and standard deviation obtained:  ", expectedValue,", ",Sqrt[variance]];
kmax=10; corr=ConstantArray[1.0,kmax+1];
Do[
corr[[k+1]]={k,(points[[1;;n-k]].points[[k+1;;n]]/(n-k)-expectedValue^2)/variance},
{k,0,kmax}
];
ListPlot[corr,PlotRange->All]

Applied to the distribution with several maxima (the choice of δ plays an important role in whether we get reasonable results):

n=100000;
a=0;b=6π;
ρ:=2((Sin[#])^2)/b+0.00001&;
x0=(a+b)/2; step=(b-a)/10;
points=myRandomMetropolis1D[n,ρ,a,b,x0,step];
ListPlot[{points},Joined->{True}]
Show[Histogram[points,{a,b,1/8},"ProbabilityDensity"],Plot[ρ[x],{x,a,b},PlotStyle->Red]]

Metropolis algorithm for n-dimensional case:

myRandomMetropolisND[n_,dim_,ρ_,a_,b_,x0_,δ_]:=Module[{x=x0,xtrial,ratio,ρx,out=ConstantArray[0.0,{n,dim}],accepted=0},
ρx=ρ[x];
Do[
(* try a new random step *)
xtrial=x+δ RandomReal[{-1,1},dim];
(* check if we are still in the given region, if not, adjust *)
Do[
If[xtrial[[d]]<a[[d]],xtrial[[d]]=a[[d]]+δ[[d]] RandomReal[]];
If[xtrial[[d]]>b[[d]],xtrial[[d]]=b[[d]]-δ[[d]] RandomReal[]],
{d,1,dim}
];
(* decide whether to make this step according to the distribution *)
ratio=ρ[xtrial]/ρx;
If[ratio⩾1||ratio>RandomReal[],x=xtrial;ρx=ρx * ratio,accepted++];
out[[i]]=x,
{i,1,n}
];
Print["Number of accepted trials: ",accepted, " from ", n];
Return[out]
];

n=10000;
a={0,0};b={2π,3π};
(* ρ:=2((Sin[#1])^2)((Sin[#2])^2)/Times@@b &; *)
ρ:=2((Sin[#[[1]]])^2)((Sin[#[[2]]])^2)/(Times@@b)+0.00001 &;
x0={3,4};δ={5,5};
points=myRandomMetropolisND[n,2,ρ,a,b,x0,δ];
ListPlot[{points},Joined->{True},PlotRange->All]
Histogram3D[points,{{a[[1]],b[[1]],1/8},{a[[2]],b[[2]],1/8}}]

Plot3D[ρ[{x,y}],{x,a[[1]],b[[1]]},{y,a[[2]],b[[2]]},PlotRange->All]

Barker algorithm

One dimensional case, ρ must be a function of x:

myRandomBarker1D[n_,ρ_,a_,b_,x0_,δ_]:=
Module[{x=x0,xtrial,ratio,ρx,ρtrial,out=ConstantArray[0.0,n],accepted=0},
ρx=ρ[x];
Do[
(* try a new random step *)
xtrial=x+δ RandomReal[{-1,1}];
(* check if we are still in the given region, if not, adjust *)
If[xtrial<a,xtrial=a+δ RandomReal[]];
If[xtrial>b,xtrial=b-δ RandomReal[]];
ρtrial=ρ[xtrial];
ratio=ρtrial/(ρx+ρtrial);
(* decide whether to make this step according to the distribution *)
If[ratio>RandomReal[],x=xtrial;ρx=ρtrial;accepted++];
out[[i]]=x,
{i,1,n}
];
Print["Number of accepted trials: ",accepted, " from ", n];
Return[out]
];

Applied to the normal distribution:

n=10^5;
a=0;b=20;
μ=10; σ=3;
ρ:=PDF[NormalDistribution[μ,σ],#]&;
x0=10.0; step=5.0;
points=myRandomBarker1D[n,ρ,a,b,x0,step];
data=RandomReal[NormalDistribution[μ,σ],{n}];
ListPlot[{data,points},Joined->{False,True}]
Show[Histogram[points,{a,b,1/4},"ProbabilityDensity"],Plot[PDF[NormalDistribution[μ,σ],x],{x,a,b},PlotStyle->Red]]
expectedValue=Total[points]/n ;
variance=Total[points^2]/n-expectedValue^2 ;
Print["Expected value and standard deviation should be: ", N[μ],", ",N[σ]];
Print["Expected value and standard deviation obtained:  ", expectedValue,", ",Sqrt[variance]];
kmax=10; corr=ConstantArray[1.0,kmax+1];
Do[
corr[[k+1]]={k,(points[[1;;n-k]].points[[k+1;;n]]/(n-k)-expectedValue^2)/variance},
{k,0,kmax}
];
ListPlot[corr,PlotRange->All]

Applied to the distribution with several maxima (the choice of δ plays an important role in whether we get reasonable results):

n=100000;
a=0;b=6π;
ρ:=2((Sin[#])^2)/b+0.00001&;
x0=(a+b)/2; step=(b-a)/10;
points=myRandomBarker1D[n,ρ,a,b,x0,step];
ListPlot[{points},Joined->{True}]
Show[Histogram[points,{a,b,1/8},"ProbabilityDensity"],Plot[ρ[x],{x,a,b},PlotStyle->Red]]

Barker algorithm for n-dimensional case:

myRandomBarkerND[n_,dim_,ρ_,a_,b_,x0_,δ_]:=Module[{x=x0,xtrial,ratio,ρx,ρtrial,out=ConstantArray[0.0,{n,dim}],accepted=0},
ρx=ρ[x];
Do[
(* try a new random step *)
xtrial=x+δ RandomReal[{-1,1},dim];
(* check if we are still in the given region, if not, adjust *)
Do[
If[xtrial[[d]]<a[[d]],xtrial[[d]]=a[[d]]+δ[[d]] RandomReal[]];
If[xtrial[[d]]>b[[d]],xtrial[[d]]=b[[d]]-δ[[d]] RandomReal[]],
{d,1,dim}
];
(* decide whether to make this step according to the distribution *)
ρtrial=ρ[xtrial];
ratio=ρtrial/(ρx+ρtrial);
If[ratio>RandomReal[],x=xtrial;ρx=ρtrial,accepted++];
out[[i]]=x,
{i,1,n}
];
Print["Number of accepted trials: ",accepted, " from ", n];
Return[out]
];

n=10000;
a={0,0};b={2π,3π};
(* ρ:=2((Sin[#1])^2)((Sin[#2])^2)/Times@@b &; *)
ρ:=2((Sin[#[[1]]])^2)((Sin[#[[2]]])^2)/(Times@@b)+0.00001 &;
x0={3,4};δ={21,2};
points=myRandomBarkerND[n,2,ρ,a,b,x0,δ];
ListPlot[{points},Joined->{True},PlotRange->All]
Histogram3D[points,{{a[[1]],b[[1]],1/8},{a[[2]],b[[2]],1/8}}]