Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Oct 2016
Lecture topic: Monte Carlo Integration & Pseudo and Quasi-Random Sequences

Integration As Area Under Curve

Here is a known function:

Clear[g]
g[x_]=Exp[x]/(1+x^2)

plot1=Plot[g[x],{x,-3,3},PlotStyle->{Red,Thick},BaseStyle->Large,ImageSize->{500,Automatic},AxesLabel->{"x","y"},AspectRatio->1/2]

The integral of this function for x∈(-1,2) is given by the area under the curve.

areaplot=Plot[g[x],{x,-1,2},PlotStyle->{Red,Thickness[0.01]},BaseStyle->Large,ImageSize->{500,Automatic},AxesLabel->{"x","y"},PlotRange->{{-3,3},{0,2}},AspectRatio->(1/2),Filling->Axis]

Show[plot1,areaplot,PlotRange->{{-3,3},{0,2}},AspectRatio->1/2]

Let's have Mathematica do the integral so we can compare to something as we go along.

integral=SetPrecision[(∫_(-1)^2)g[x]\[DifferentialD]x,12]//Re

Integration as Ratio of Points

The area can be thought of as a certain percentage of the area given in the rectangle, where the width of the rectangle is chosen as the integration limits and the height contains the entire area that represents the integral. For this particular curve, I have chosen the height to be g(2).

rectangle=Plot[{g[2],0},{x,-1,2},ImageSize->{500,Automatic},AxesLabel->{"x","y"},PlotRange->{{-3,3},{0,2}},AspectRatio->(1/2),Filling->{1->{2}}];

mainplot=Show[rectangle,plot1,areaplot,Axes->True,AspectRatio->1/2]

Let’s set up a grid of points, and count how many points are below the curve.

point[x_,y_]:=Graphics[{Blue,Disk[{x,y},{0.05,0.05*2*2/6}]}]
listxy={};
Do[listxy=Append[listxy,point[x,y]],{x,-1,2,0.3},{y,0,g[2],0.3}]
Show[mainplot,listxy,AspectRatio->1/2,PlotRange->{All,All}]

Therefore,

(∫_(-1)^2)g(x)dx = area under curve
~((number of points below g(x))/(number of points used))×area of blue rectangle
~(44/55)×(3*g(2))

SetPrecision[(44/55)*(3*g[2]),12]

Compare to the actual result:

integral

Not so good. However, we only used 55 points. Let's use more!

listxy={};
Do[listxy=Append[listxy,point[x,y]],{x,-1,2,0.09},{y,0,g[2],0.09}];
Show[mainplot,listxy,AspectRatio->1/2,PlotRange->{All,All}]

The graphics are not helping us anymore, so we will drop them. The counting process can be automated. Now, rather than set the number of points we will set a grid, based on widths Δx and Δy.

g[x_]=Exp[x]/(1+x^2);
a=-1;
b=2;
c=0;
d=g[2];
area=(b-a)*(d-c);

Λx=0.1;
Λy=0.1;

pointsbelow=0;
totalpoints=0;
Do[
totalpoints=totalpoints+1;
If[g[x]>y,pointsbelow=pointsbelow+1];
,{x,a,b,Λx},{y,c,d,Λy}];

SetPrecision[(pointsbelow/totalpoints)*area,12]

totalpoints

Better. But what if we wanted something more accurate? We would need to do the entire calculation over again, with a new grid of points! That can be time consuming, but can be done.

g[x_]=Exp[x]/(1+x^2);
a=-1;
b=2;
c=0;
d=g[2];
area=(b-a)*(d-c);

Λx=0.01;
Λy=0.01;

pointsbelow=0;
totalpoints=0;
Do[
totalpoints=totalpoints+1;
If[g[x]>y,pointsbelow=pointsbelow+1];
,{x,a,b,Λx},{y,c,d,Λy}];

SetPrecision[(pointsbelow/totalpoints)*area,12]

totalpoints

The error in using 44548 points with the grid method is

SetPrecision[Abs[integral-(pointsbelow/totalpoints)*area],12]

There are two significant problems with the grid method:

To get a more accurate estimate, you need to redo the entire calulation using a new grid. This is inefficient.

If you extend this to higher dimensions, it requires a great many nested loops.

Pseudo-Random Numbers & Monte Carlo Integration

Random numbers are fairly easy to get a hold of, although the details of the techniques that leads to them can be a bit daunting. Most people simply use their favourite random number generator and don't worry too much about it. This is usually OK, unless you are using a tremendous number of points from one run of your random number generator.

Here is Mathematica's random number generator. The SeedRandom simply means the random number sequence will always start at the same point.

SeedRandom[2]
RandomReal[{-2,2}]

With random numbers, we can do the calculations above in a manner that extends easily to higher dimension, and allows us to add more points easily (which were the two problems with using the grid). Let's see how the previous problem looks when we use random numbers instead of a grid.

g[x_]=Exp[x]/(1+x^2);

(*random numbers to generate the points*)
xpoint:=RandomReal[{-1,2}]
ypoint:=RandomReal[{0,g[2]}]

(*create a graphic blue dot we can use to see where the point is*)
xypoint:=Graphics[{Blue,Disk[{xpoint,ypoint},{0.05,0.05*2*2/6}]}]

(*create a list of 100 random blue dots*)
listxy={};
Do[listxy=Append[listxy,xypoint],{numpoints,0,100}];

(*plot the situation*)
Show[mainplot,listxy,AspectRatio->1/2]

Notice there is some clustering of the points--that is a normal feature of random numbers! Now, let's do the integral:

xpoint:=RandomReal[{-1,2}]
ypoint:=RandomReal[{0,g[2]}]
g[x_]=((E^x)/(1+(x^2)));
a=-1;
b=2;
c=0;
d=g[2];
area=(b-a) (d-c);
pointsbelow=0;
totalpoints=5000;

Do[
If[g[xpoint]>ypoint,pointsbelow=pointsbelow+1];
,{ipoints,1,totalpoints}];

SetPrecision[((pointsbelow area)/totalpoints),12]
totalpoints

integral

Whoa--not baaaad!! Let's look at the error in the technique. Now, since we are using random numbers instead of a grid of points, we can get an estimate for the integral after each individual point is used. Contrast this with the grid method, where you would have to redo all the points in the grid!

xpoint:=RandomReal[{-1,2}]
ypoint:=RandomReal[{0,g[2]}]
g[x_]=((E^x)/(1+(x^2)));
a=-1;
b=2;
c=0;
d=g[2];
area=(b-a) (d-c);
pointsbelow=0;
currentpoints=0;
totalpoints=5000;

approx={};
Do[
currentpoints=currentpoints+1;
If[g[xpoint]>ypoint,pointsbelow=pointsbelow+1];approx=Append[approx,{currentpoints,SetPrecision[Abs[integral-((pointsbelow*area)/currentpoints)],12]}];,{ipoints,1,totalpoints}];

ListPlot[approx,PlotStyle->{Red},BaseStyle->Large,ImageSize->{500,Automatic},Frame->True,FrameLabel->{"Number of Points","Error"},PlotRange->{All,{0,0.1}}]

Let's redo the calculation with more points:

xpoint:=RandomReal[{-1,2}]
ypoint:=RandomReal[{0,g[2]}]
g[x_]=((E^x)/(1+(x^2)));
a=-1;
b=2;
c=0;
d=g[2];
area=(b-a) (d-c);
pointsbelow=0;
currentpoints=0;
totalpoints=15000;

approx={};
Do[
currentpoints=currentpoints+1;
If[g[xpoint]>ypoint,pointsbelow=pointsbelow+1];approx=Append[approx,{currentpoints,SetPrecision[Abs[integral-((pointsbelow*area)/currentpoints)],12]}];,{ipoints,1,totalpoints}];

ListPlot[approx,PlotStyle->{Red},BaseStyle->Large,ImageSize->{600,Automatic},Frame->True,FrameLabel->{"Number of Points","Error"},PlotRange->{All,{0,0.1}}]

Note that the error goes roughly as (n^(-1/2)), where n is the number of points used. Of course, since we are using random numbers, the actual result does not follow the predicted result exactly. Let’s plot the error on a logarithmic scale.

errorplot2=LogLogPlot[(n^(-1/2)),{n,1,15000},PlotStyle->{Blue,Thick}];
errorplot=ListLogLogPlot[approx,PlotStyle->{Red},Frame->True,ImageSize->{600,Automatic}];

Show[errorplot2,errorplot,AxesLabel->{"Number of Points","Error"},Frame->True,BaseStyle->20]

Corput Sequence to Generate quasi-random numbers

Pseudo vs Quasi Random Numbers

Here are the random numbers Mathematica is generating and how they fill space in (ℝ^2) :

SeedRandom[2];
random=Table[RandomReal[{0,1},2],{i,1,2000}];

rplot=ListPlot[random,Frame->True,AspectRatio->1,PlotStyle->Blue,ImageSize->{500,Automatic}]

Again, notice that there is some clustering, which we expect with random numbers. This means there are some places that are not sampled at all, which is bad if something interesting is happening in that region. In fact, Mathematica is using pseudo-random numbers (deterministically generated numbers that strive to be random in the sense of being independent and identically distributed).

For numerical integration (which is what we are focussed on) psuedo-random numbers are not the best choice to use. We don’t want a grid of numbers for the reasons mentioned earlier, but using psuedo-random numbers give estimates to integrals that converge slowly. What we really want is to efficiently fill space in some manner that is not grid-like, but is more regular than pseudo-random where each point is independent. We can achieve this if we exchange the requirement of independence with self-avoiding.

This leads us to look for quasi-random sequences - which really aren't random at all! They use some mathematical technique to generate a sequence of points that fill phase space but avoid clustering. The sobol sequence is possibly the most useful quasi-random sequence. Here, we construct another, simpler quasi-random sequence called the Corput sequence.

Quasi-Random: Corput Sequence in (ℝ^1)

Recall the base expansion for a number:

928.64=9×(10^2)+2×(10^1)+8×(10^0)+6×(10^(-1))+4×(10^(-2))

Here is a Mathematica commands you may not have seen before, which help us pull out the coefficients:

IntegerDigits[928,10]

This works in any base:

IntegerDigits[6,2]

The Corput sequence is defined as follows:

The decimal number n is expanded in base b.

The number in base b is reflected.

Convert back to decimal number.

where we do this process for n=1,2,3,4,... to generate the Corput sequence. Here is the Corput number calculated when n=6 and b=2:

Expand in base 2: 6=1×(2^2)+1×(2^1)+0×(2^0)-> 110 in base 2
Reflect in base 2: 011-> 0×(2^(-1))+1×(2^(-2))+1×(2^(-3))
Convert to decimal: (1/2)+(1/8)=(3/8)

This module automates the process outlined by hand:

Corput[n_,b_]:=Module[{expandbaseb,m},
expandbaseb=IntegerDigits[n,b];
m=Length[expandbaseb];
expandbaseb=Reverse[expandbaseb];
Return[(∑_(j=0)^(m-1))expandbaseb[[j+1]]*(b^(-j-1))];
]

Corput[6,2]

And here is the Corput sequence when the base is b=2.

corput=Table[Corput[i,2],{i,1,20}]

Here is a plot of how the Corput sequence with base 2 fills the region x ϵ (0,1). Notice that any “holes” in the region are filled by the sequence, which is what makes it a quasi-random sequence (it avoids clumping).

Clear[list]
n=80;

list={};
Do[
list=Append[list,Graphics[Disk[{Corput[i,2],0},{0.01,0.01}]]]
,{i,1,n}];

Manipulate[
Show[Take[list,i],PlotRange->{{0,1},{-(1/8),(1/8)}},Axes->{True,False},ImageSize->{500,Automatic},AspectRatio->(1/4)]
,{i,1,n,1}]

Here is a plot of how the Corput sequence with base 3 fills the region x ϵ (0,1).

corput=Table[Corput[i,3],{i,1,20}]

Clear[list2]
n=80;

list2={};
Do[
list2=Append[list2,Graphics[Disk[{Corput[i,3],0},{0.01,0.01}]]]
,{i,1,n}];

Manipulate[
Show[Take[list2,i],PlotRange->{{0,1},{-(1/8),(1/8)}},Axes->{True,False},ImageSize->{500,Automatic},AspectRatio->(1/4)]
,{i,1,n,1}]

Quasi-Random: Halton Sequence in (ℝ^2)

To extend to higher dimensions, you take as base the next prime number, so in (ℝ^2) you use the x-coordinate with base 2, the y-coordinate with base 3.

Halton2D[n_]:={Corput[n,2],Corput[n,3]}

halton=Table[Halton2D[k],{k,1,10}]

Here are some examples of how the Halton sequence fills space in (ℝ^2).

halton=Table[Halton2D[k],{k,1,2000}];

hplot=ListPlot[halton,Frame->True,AspectRatio->1,PlotStyle->Red,ImageSize->{500,Automatic}]

Here is an animation of how the Halton sequence fills the x y-plane.

Clear[list3]
n=1200;

list3={};
Do[
list3=Append[list3,Graphics[{Red,Disk[Halton2D[i],{0.01,0.01}]}]]
,{i,1,n}];

Manipulate[
Show[Take[list3,i],PlotRange->{{0,1},{0,1}},Frame->True,AspectRatio->1,ImageSize->{500,Automatic}]
,{i,1,n,1}]

Quasi-Random: Halton Sequence in (ℝ^3)

To extend to higher dimensions, you take as base the next prime number, so in (ℝ^3) you use the x-coordinate with base 2, the y-coordinate with base 3 and the z-coordinate with base 5. This is called a Halton sequence.

Halton3D[n_]:={Corput[n,2],Corput[n,3],Corput[n,5]}

halton=Table[Halton3D[k],{k,1,10}]

Here are some examples of how the Halton sequence fills 3D space. Rotate the image and notice that there is no clumping, the hallmark of a good quasi-random sequence.

halton=Table[Halton3D[k],{k,1,2000}];

hplot3D=ListPointPlot3D[halton,ImageSize->{500,Automatic},PlotStyle->Purple]

May as well animate!

Clear[list4]
n=200;

list4={};
Do[
list4=Append[list4,Graphics3D[{Purple,Sphere[Halton3D[i],0.01]}]]
,{i,1,n}];

Manipulate[
Show[Take[list4,i],PlotRange->{{0,1},{0,1},{0,1}},ImageSize->{500,Automatic}]
,{i,1,n,1}]

Comparison of Pseudo-Random to Quasi-Random

Show[rplot,hplot,ImageSize->{600,Automatic}]

Notice that the Halton sequence does not have the clumping that a true random (or pseudo-random) sequence possesses. This makes the quasi-random Halton sequence useful for sampling space, and doing integrals!

Mathematica’s Built-In Random Number Generators

(ℝ^2)

You can access a variety of different psuedo-random generators using the SeedRandom command in Mathematica. Two of the ones below are quasi-random, not pseudo-random. Can you guess which two are quasi-random?

ListPlot[BlockRandom[SeedRandom[1];
RandomReal[1,{2000,2}]],AspectRatio->1]

ListPlot[BlockRandom[SeedRandom[1,Method->"Rule30CA"];
RandomReal[1,{2000,2}]],AspectRatio->1]

ListPlot[BlockRandom[SeedRandom[1,Method->"MersenneTwister"];
RandomReal[1,{2000,2}]],AspectRatio->1]

ListPlot[BlockRandom[SeedRandom[Method->{"MKL",Method->{"Niederreiter","Dimension"->2}}];
RandomReal[1,{2000,2}]],AspectRatio->1]

ListPlot[BlockRandom[SeedRandom[Method->{"MKL",Method->{"Sobol","Dimension"->2}}];
RandomReal[1,{2000,2}]],AspectRatio->1]

ListPlot[BlockRandom[SeedRandom[1,Method->{"ParallelMersenneTwister","Index"->2}];
RandomReal[1,{2000,2}]],AspectRatio->1]

(ℝ^3)

ListPointPlot3D[BlockRandom[SeedRandom[1];RandomReal[1,{1000,3}]],AspectRatio->1,PlotStyle->{Purple}]

ListPointPlot3D[BlockRandom[SeedRandom[1,Method->"Rule30CA"];RandomReal[1,{1000,3}]],AspectRatio->1,PlotStyle->{Purple}]

ListPointPlot3D[BlockRandom[SeedRandom[1,Method->"MersenneTwister"];RandomReal[1,{1000,3}]],AspectRatio->1,PlotStyle->{Purple}]

ListPointPlot3D[BlockRandom[SeedRandom[Method->{"MKL",Method->{"Niederreiter","Dimension"->3}}];RandomReal[1,{1000,3}]],AspectRatio->1,PlotStyle->{Purple}]

ListPointPlot3D[BlockRandom[SeedRandom[Method->{"MKL",Method->{"Sobol","Dimension"->3}}];RandomReal[1,{1000,3}]],AspectRatio->1,PlotStyle->{Purple}]

ListPointPlot3D[BlockRandom[SeedRandom[1,Method->{"ParallelMersenneTwister","Index"->2}];
RandomReal[1,{1000,3}]],AspectRatio->1,PlotStyle->{Purple}]

Extend Integration to Double Integral

A Double Integral

In higher dimension, the technique can be used to represent a volume (double integral), which is really were Monte Carlo integration starts to shine.

(∫_(-1)^2)(∫_(-1)^2)g(x,y)dx dy
~((number of points below g(x,y))/(number of points used))×volume points are samples from

g[x_,y_]=((Exp[x])/((1+(x^2))))((Exp[y])/((2+(y^2))))

Plot3D[g[x,y],{x,-1,2},{y,-1,2}]

Here is the result for us to compare to:

integral=NIntegrate[g[x,y],{x,-1,2},{y,-1,2}]

Mathematica’s Default Pseudo-Random Numbers

Here we set up an approximation to the integral using Mathematica’s psuedo-random numbers. Since we are using a lot of points, we are only storing every 100th point for plotting purposes.

Clear[xpoint,ypoint,zpoint]
SeedRandom[3]
x1=-1;x2=2;
y1=-1;y2=2;
z1=0;z2=g[2,2];
volume=(x2-x1) (y2-y1) (z2-z1);
xpoint:=RandomReal[{x1,x2}]
ypoint:=RandomReal[{y1,y2}]
zpoint:=RandomReal[{z1,z2}]
pointsbelow=0;
currentpoints=0;
totalpoints=50000;
approx={};
Do[
currentpoints=currentpoints+1;
If[g[xpoint,ypoint]>zpoint,pointsbelow=pointsbelow+1];If[Floor[(ipoints/100)] 100==ipoints,AppendTo[approx,{currentpoints,SetPrecision[Abs[integral-((pointsbelow volume)/currentpoints)],12]}]];,{ipoints,1,totalpoints}];

errorplotrandom=ListLogLogPlot[approx,PlotStyle->{Red},ImageSize->{500,Automatic},AxesLabel->{"Number of Points","Error"},PlotRange->{All,All},BaseStyle->Medium]

Our Corput Quasi-Random Numbers

Here we set up an approximation to the integral using a Halton sequence. Notice that since the original Halton sequence was on the interval [0,1] we are transforming to a laerger integral based on the limits of integration. Again, since we are using a lot of points, we are only storing every 100th point for plotting purposes.

Corput[n_,b_]:=Module[{expandbaseb,m},
expandbaseb=IntegerDigits[n,b];
m=Length[expandbaseb];
expandbaseb=Reverse[expandbaseb];
Return[(∑_(j=0)^(m-1))expandbaseb[[j+1]]*(b^(-j-1))];
]

Clear[xpoint,ypoint,zpoint]

x1=-1;x2=2;
y1=-1;y2=2;
z1=0;z2=g[2,2];
volume=(x2-x1)*(y2-y1)*(z2-z1);

xpoint[n_]:=x1+(x2-x1)*Corput[n,2]//N
ypoint[n_]:=y1+(y2-y1)*Corput[n,3]//N
zpoint[n_]:=z1+(z2-z1)*Corput[n,5]//N

pointsbelow=0;
currentpoints=0;
totalpoints=50000;

approx={};
Do[
currentpoints=currentpoints+1;
If[g[xpoint[ipoints],ypoint[ipoints]]>zpoint[ipoints],pointsbelow=pointsbelow+1];
If[Floor[ipoints/100]*100==ipoints,approx=Append[approx,{currentpoints,SetPrecision[Abs[integral-(pointsbelow/currentpoints)*volume],12]}]];
,{ipoints,1,totalpoints}];

errorplotquasirandom=ListLogLogPlot[approx,PlotStyle->{Blue},ImageSize->{500,Automatic},AxesLabel->{"Number of Points","Error"},PlotRange->{All,All},BaseStyle->Medium]

Mathematica’s Sobol Quasi-Random Numbers

Here we set up an approximation to the integral using Mathematica’s implementation of sobol quasi-random numbers. Since we are using a lot of points, we are only storing every 100th point for plotting purposes.

Clear[xpoint,ypoint,zpoint]
SeedRandom[Method->{"MKL",Method->{"Sobol","Dimension"->3}}]
x1=-1;x2=2;
y1=-1;y2=2;
z1=0;z2=g[2,2];
volume=(x2-x1) (y2-y1) (z2-z1);
xpoint:=RandomReal[{x1,x2}]
ypoint:=RandomReal[{y1,y2}]
zpoint:=RandomReal[{z1,z2}]
pointsbelow=0;
currentpoints=0;
totalpoints=50000;
approx={};
Do[
currentpoints=currentpoints+1;
If[g[xpoint,ypoint]>zpoint,pointsbelow=pointsbelow+1];If[Floor[(ipoints/100)] 100==ipoints,approx=Append[approx,{currentpoints,SetPrecision[Abs[integral-((pointsbelow volume)/currentpoints)],12]}]];,{ipoints,1,totalpoints}];

errorplotsobol=ListLogLogPlot[approx,PlotStyle->{Darker[Green]},ImageSize->{500,Automatic},AxesLabel->{"Number of Points","Error"},PlotRange->{All,All},BaseStyle->Medium]

Comparison

Let’s compare the errors:

Show[errorplotrandom,errorplotquasirandom,errorplotsobol]

Notice that the quasi-random sequences provide a better estimate for smaller n values (although in this example, by the time we get out to n=40000 they are both about the same). The text discusses the error in quasi-random, which goes as (n^(-1/2-1/(2d))) or (((log(n)))^d)(n^(-1)) where d is the dimension of the random numbers

errorplot=LogLogPlot[{(nn^(-1/2)),(nn^(-1/2-1/8)),(((Log[nn]))^2)(nn^(-1))},{nn,1,50000},PlotStyle->{{Red,Thick},{Blue,Thick},{Green,Thick}}];
Show[errorplot,errorplotrandom,errorplotquasirandom,errorplotsobol]

The error is defined in terms of an average, so for any single computation we would expect the result to vary. We can get a better understanding of the error if we instead fix n and run a series of trials and take the average error over the trials. This is what is done in the text.

A Bit More In Depth Exploration of the Error

The error is defined in terms of an average, so for any single computation we would expect the result to vary. We can get a better understanding of the error if we instead fix n and run a series of trials and take the average error over the trials.

First, do for a pseudo-random sequence:

Clear[xpoint,ypoint,zpoint]

(*definte the integral we are investigating*)
g[x_,y_]=((Exp[x])/((1+(x^2))))((Exp[y])/((2+(y^2))));
x1=-1.0;
x2=2.0;
y1=-1.0;
y2=2.0;
z1=0.0;
z2=g[2.0,2.0];
volume=(x2-x1) (y2-y1) (z2-z1);
integral=Integrate[g[x,y],{x,x1,x2},{y,y1,y2}];

(*do for different totalpoints*)
errormeanpseudo={};
Do[
(*do 10 trials at a given number of totalpoints*)
error={};
Do[
SeedRandom[numrun*totalpoints];
pointsbelow=0;
currentpoints=0;

(*compute the approximation to the integral for given trial*)
Do[
xpoint=RandomReal[{x1,x2}];
ypoint=RandomReal[{y1,y2}];
zpoint=RandomReal[{z1,z2}];
currentpoints=currentpoints+1;
If[g[xpoint,ypoint]>zpoint,pointsbelow=pointsbelow+1];
,{ipoints,1,totalpoints}];
(*store the error for given trial*)
error=Append[error,SetPrecision[Abs[integral-((pointsbelow*volume)/currentpoints)],12]];
,{numrun,1,10,1}];

(*store the mean error for given totalpoints*)
Print[totalpoints," ",Mean[error]];
errormeanpseudo=Append[errormeanpseudo,{totalpoints,Mean[error]}];
,{totalpoints,{100,200,500,1000,2000,3000,5000,7000,10000,15000,20000}}]

Now do for a quasi-random sequence (make sure the Corput is defined!). This is a bit slower than the previous one.

Corput[n_,b_]:=Module[{expandbaseb,m},
expandbaseb=IntegerDigits[n,b];
m=Length[expandbaseb];
expandbaseb=Reverse[expandbaseb];
Return[(∑_(j=0)^(m-1))expandbaseb[[j+1]]*(b^(-j-1))];
]

Clear[xpoint,ypoint,zpoint]

(*definte the integral we are investigating*)
g[x_,y_]=((Exp[x])/((1+(x^2))))((Exp[y])/((2+(y^2))));
x1=-1.0;
x2=2.0;
y1=-1.0;
y2=2.0;
z1=0.0;
z2=g[2.0,2.0];
volume=(x2-x1) (y2-y1) (z2-z1);
integral=Integrate[g[x,y],{x,x1,x2},{y,y1,y2}];

(*do for different totalpoints*)
errormeanquasi={};
Do[
(*do 10 trials at a given number of totalpoints*)
error={};
Do[
SeedRandom[numrun*totalpoints];
pointsbelow=0;
currentpoints=0;

(*compute the approximation to the integral for given trial*)
Do[
xpoint=x1+(x2-x1)*Corput[ipoints,2]//N;
ypoint=y1+(y2-y1)*Corput[ipoints,3]//N;
zpoint=z1+(z2-z1)*Corput[ipoints,5]//N;
currentpoints=currentpoints+1;
If[g[xpoint,ypoint]>zpoint,pointsbelow=pointsbelow+1];
,{ipoints,1,totalpoints}];
(*store the error for given trial*)
error=Append[error,SetPrecision[Abs[integral-((pointsbelow*volume)/currentpoints)],12]];
,{numrun,1,10,1}];

(*store the mean error for given totalpoints*)
Print[totalpoints," ",Mean[error]];
errormeanquasi=Append[errormeanquasi,{totalpoints,Mean[error]}];
,{totalpoints,{100,200,500,1000,2000,3000,5000,7000,10000,15000,20000}}]

errorplots=ListLogLogPlot[{errormeanpseudo,errormeanquasi},PlotStyle->{{Red,PointSize[0.02]},{Blue,PointSize[0.02]}},ImageSize->{500,Automatic},AxesLabel->{"Number of Points","Error"},PlotRange->{All,All},BaseStyle->Large]

The quasi-random is better for integration, as it will generally give better results for smaller n that pseudo-random.

Potential Problems in Higher Dimension

Quasi-Random Failure: Halton Sequence in (ℝ^232)

I've said Corput doesn't work well for high dimensional random numbers. Let check that out. I will work with decimal to speed things up, but you would see similar behaviour if you worked with integers (the behaviour we see is not due to round-off error).

Corput[n_,b_]:=Module[{expandbaseb,m},
expandbaseb=IntegerDigits[n,b];
m=Length[expandbaseb];
expandbaseb=Reverse[expandbaseb];
Return[N[(∑_(j=0)^(m-1))expandbaseb[[j+1]]*(b^(-j-1))]];
]

Let’s get a table of the first 232 primes, which will be used for each dimension in the Halton sequence.

dim=232;
base=Table[Prime[i],{i,1,dim}]

This defines the Halton sequence in (ℝ^232) that we can investigate. I’ve output the first vector in the sequence.

Halton[n_]:=Table[Corput[n,base[[i]]],{i,1,dim}]
Halton[1]

This creates a set of 4000 points in (ℝ^232) that we can investigate. The quantity “halton” contains 4000 vectors of length 232. This will take a little while to create.

npts=4000;
halton=Table[Halton[k],{k,1,npts}];

Length[halton]
Length[halton[[1]]]

Now, let’s choose to look at a cross section of the (ℝ^232) space.

points=Table[{halton[[i,2]],halton[[i,193]]},{i,1,npts}];

Now, we know that quasi-random numbers should avoid clumping, so the sketch in (ℝ^2) should look similar to what we saw earlier.Notice that we have visual evidence that the halton sequence does NOT work well in high dimension--look at the bands that are appearing in our cut of dimensions 2 and 193. It is not sampling this high dimension space in the way we want for a quasi-random sequence.

hplot=ListPlot[points,Frame->True,AspectRatio->1,PlotStyle->Red,ImageSize->{500,Automatic}]

The visual evidence is even more striking if you examine a slice in (ℝ^3) (compare to Figure 9.3 in the text). Rotate the plot below to be amazed and saddened! The Halton sequence is worthless for 232 dimensional random numbers!

points3D=Table[{halton[[i,2]],halton[[i,163]],halton[[i,200]]},{i,1,npts}];

hplot3D=ListPointPlot3D[points3D,ImageSize->{500,Automatic},PlotStyle->Purple]

If you need to use a high dimension random number generator, you should step up to a sobol sequence which is good to 32,000 dimensions. Halton is most likely safe for under 14 dimensions. Of course, if you are doing anything mission critical just use a sobol sequence (we will see how well it does integrals shortly)!

Pseudo-Random Success in (ℝ^232)

We’ve seen that Corput doesn't work well for high dimensional random numbers. Let check sobol out, which is supposedly better for higher dimensional randum vectors.

Now, let’s choose to look at a cross section of the (ℝ^232) space. First, let’s see how Mathematica can be used to create a set of random vectors that is small so we can see the output! Note that I had trouble making this work until I choose to initialize with a different method. Not sure why that was the case.

dim=12;
npts=5;
SeedRandom[1,Method->"Rule30CA"]

points={};
Do[
AppendTo[points,RandomReal[{-1,1},dim]];
,{i,1,npts}];

slice2D={};
slice3D={};
Do[
AppendTo[slice2D,{points[[i,2]],points[[i,7]]}];
AppendTo[slice3D,{points[[i,3]],points[[i,7]],points[[i,12]]}];
,{i,1,npts}];

Here are the slices we took:

points//MatrixForm
slice2D//MatrixForm
slice3D//MatrixForm

Now that we are satisfied it is working, make it huge:

dim=232;
npts=4000;
SeedRandom[1,Method->"Rule30CA"]

points={};
Do[
points=Append[points,RandomReal[{0,1},dim]];
,{i,1,npts}];

slice2D={};
slice3D={};
Do[
slice2D=Append[slice2D,{points[[i,2]],points[[i,7]]}];slice3D=Append[slice3D,{points[[i,3]],points[[i,12]],points[[i,50]]}];
,{i,1,npts}];

Here are the slices, and note that each seems to have a decent spread of points.

ListPlot[slice2D]

ListPointPlot3D[slice3D,ImageSize->{500,Automatic},PlotStyle->Purple]

This looks good in terms of how space is being sampled.

Pros/Cons of Monte Carlo and Quasi-Monte Carlo Integration

The main benefits of Monte-Carlo integration are for integrals of three types:

Integrals over a complicated domain,

Multidimensional integrals (which are difficult by any means!), although you have to be careful when using quasi-monte carlo methods and dimensions over about 16 (the lack of randomness in higher dimensions causes some problems with the technique),

Integrals over integrands which are essentially constant, or essentially Gaussian in nature.

The main problem with Monte Carlo integration is that it can take a huge number of points to get convergence. Therefore, it is best used when high accuracy is not desired. For these reasons, it is used frequently in semiclassical quantum chemistry calculations, or other molecular dynamics calculations.

What proportion of symmetric 2×2 matrices with entries in [0,1] have positive determinant? A Monte Carlo Solution Two Ways

Solution 1

Let’s first just pick points in C:[0,1]×[0,1]×[0,1] and compute if x y-(z^2)>0 where our generic matrix is (({{x, z}, {z, y}})).

Corput[n_,b_]:=Module[{expandbaseb,m},
expandbaseb=IntegerDigits[n,b];
m=Length[expandbaseb];
expandbaseb=Reverse[expandbaseb];
Return[(∑_(j=0)^(m-1))expandbaseb[[j+1]]*(b^(-j-1))];
]

We just modify the Monte Carlo process to count the number of points for which x y-(z^2)>0.

Clear[xpoint,ypoint,zpoint]

x1=0;
x2=1;
y1=0;
y2=1;
z1=0;
z2=1;

pointscount=0;
currentpoints=0;
totalpoints=100000;

approx={};
Do[
xpoint=x1+(x2-x1)*Corput[ipoints,2]//N;
ypoint=y1+(y2-y1)*Corput[ipoints,3]//N;
zpoint=z1+(z2-z1)*Corput[ipoints,5]//N;
currentpoints=currentpoints+1;

If[xpoint*ypoint-(zpoint^2)>0,pointscount=pointscount+1];
If[Floor[ipoints/100]*100==ipoints,approx=Append[approx,{currentpoints,SetPrecision[(pointscount/currentpoints),12]}]];
,{ipoints,1,totalpoints}];

plotquasirandom=ListLogLogPlot[approx,PlotStyle->{Blue},ImageSize->{500,Automatic},AxesLabel->{"Number of Points","Proportion"},PlotRange->{All,All}]

Here is our last approximation:

Last[approx]

Solution 2

If we want x y-(z^2)>0, we can look at an integral over the region. In fact, setting this up will show us what the exact answer must be, rendering the Monte-Carlo process moot!

x y > (z^2)
0 < (z^2)< x y
0 < z <(Sqrt[x y])

So the integral that tells us the answer is

(∫_0^1)(∫_0^1)(∫_0^(Sqrt[x y]))1dzdxdy

Integrate[Integrate[1,{z,0,Sqrt[x y]}],{x,0,1},{y,0,1}]

What we found above was close to this. Notice the Monte-Carlo routine is virtually identical to what we did above.

(∫_0^1)(∫_0^1)(Sqrt[x y])dxdy

Corput[n_,b_]:=Module[{expandbaseb,m},
expandbaseb=IntegerDigits[n,b];
m=Length[expandbaseb];
expandbaseb=Reverse[expandbaseb];
Return[(∑_(j=0)^(m-1))expandbaseb[[j+1]]*(b^(-j-1))];
]

Clear[xpoint,ypoint,zpoint]

x1=0;
x2=1;
y1=0;
y2=1;
z1=0;
z2=1;

pointscount=0;
currentpoints=0;
totalpoints=100000;

approx={};
Do[
xpoint=x1+(x2-x1)*Corput[ipoints,2]//N;
ypoint=y1+(y2-y1)*Corput[ipoints,3]//N;
zpoint=z1+(z2-z1)*Corput[ipoints,5]//N;
currentpoints=currentpoints+1;

If[zpoint<(Sqrt[xpoint*ypoint]),pointscount=pointscount+1];
If[Floor[ipoints/100]*100==ipoints,approx=Append[approx,{currentpoints,SetPrecision[(pointscount/currentpoints),12]}]];
,{ipoints,1,totalpoints}];

plotquasirandom=ListLogLogPlot[approx,PlotStyle->{Blue},ImageSize->{500,Automatic},AxesLabel->{"Number of Points","Proportion"},PlotRange->{All,All}]

Here is our last approximation:

Last[approx]

A Mathematica Module using Sobol Sequence for Roots of Unity (Revisiting!)

Each point is space can be used to initialize Newton’s method to find the root of a complex valued equation. In this example, the equation is (x^n)-((Subscript[c, 0])+(Subscript[c, 1]) i)=0 for a=(Subscript[c, 0])+(Subscript[c, 1])i complex

This module uses the sobol quasi-random sequence. It produces a much prettier graphic at the end than if we just used pseudo-random numbers, since the points do not suffer from clustering. It is also faster than our Halton sequence (this should only take a minute or two to run).

SeedRandom[2,Method->{"MKL",Method->{"Sobol","Dimension"->2}}]

This structure returns one two dimensional randum number on the interval (-1,1):

hold=RandomReal[{-1,1},{1,2}]

hold[[1,1]]

plotrootysobol[c0_,c1_,n_,npts_]:=Module[{y,r,f,fp,hold},
y=c0+c1 I;
r=(Sqrt[(c0^2)+(c1^2)]);
theta=ArcTan[(c1/c0)];
root[i_]:=N[(r^(1/n)) (Cos[((theta+2 π i)/n)]+I Sin[((theta+2 π i)/n)])];
f[x_]:=(x^n)-y;
fp[x_]:=f'[x];
Do[set[i]={},{i,0,n-1}];
acc=4;
tol=(10.^(-acc));
max=500;
Do[Clear[p0,p1,i,p,q,a,b];
hold=RandomReal[{-1,1},{1,2}];
a=hold[[1,1]];
b=hold[[1,2]];
p0=a+b I;
i=1;
While[i≤max,
p=p0-((f[p0])/(fp[p0]));
If[Abs[f[p]]<tol,{Do[If[Abs[p-root[k]]<(10^(-acc+1)),set[k]=Append[set[k],{a,b}]],{k,0,n-1}],Break[]}];p0=p;i=i+1],{jj,1,npts}];
setlist={};
Do[setlist=Append[setlist,set[i]],{i,0,n-1}] ;
huelist={};
Do[huelist=Append[huelist,Hue[(i/n)]],{i,0,n-1}];

ListPlot[setlist,AspectRatio->1,Axes->True,Frame->True,PlotRange->{{-1,1},{-1,1}},PlotStyle->huelist,ImageSize->{700,700}]]

plot2=plotrootysobol[1.0,0.0,4,60000]

Here’s another example:

plot2=plotrootysobol[4.0,-3.0,6,100000]