In a previous post, I promised to write about testing the quality of pseudo-random number sequences.  I'll post later about some of the statistical tests often used, but I first wanted to mention a sort of practical test one can do. One of the many things you might want to do with pseudorandomly generated numbers is Monte Carlo integration/simulatation/etc.  As mentioned by acer in this comment, Monte Carlo integration can be shown to work better with some of the pseudorandom number generators (PRNGs) which are considered inferior in a statistical sense.  In this post, we will play with a simple Monte Carlo approximation of π.

In my old blog post on pseudo-random numbers I was considering sequences of random bits, but we'll need larger numbers to do a Monte Carlo approximation.  I chose to use 8-bit integers paired up into points in the square [0,255] x [0,255] and formed them using the ?Bits package.

(**) r1 := rand(0..1);
(**) size := 8;
(**) r := 2^size - 1;
(**) L := [seq(r1(), i=1..size*2*1000)]: # enough bits for 1000 8-bit points
(**) pointList := [seq( [ Bits:-Join([seq(L[j], j=size*i+1..size*i+size)]),
(**)                    Bits:-Join([seq\                                   
(**) (L[j], j=size*(i+1)+1..size*(i+1)+size)]) ]
(**)           , i=0..floor(nops(L)/size)-2, 2)]:

From the list of points, we can approximate Pi by counting the points the landed inside the circle inscribed in the square (approximating the area of the circle) and dividing by the total number of points (approximating the area of the square): this ratio should be, approximately (π*(255/2)^2) / 255^2 = π/4 if the numbers were random enough.

(**) M := map(x -> evalb(((x[1]-r/2.)^2+(x[2]-r/2.)^2) <= r^2/4.), pointList):
(**) numC := numboccur(M,true):
(**) PiApprox := evalf( 4 * numC / nops(M) );
#                            3.1960

It turns out this is not a great approximation: not quite two digits of accuracy. In fact, even if you use all 65,536 8-bit points you can get only two full digits of accuracy: 3.1152. It does make for a nice plot, however:

(**) (L1,L2) := selectremove(x -> evalb(((x[1]-r/2.)^2+(x[2]-r/2.)^2) <= r^2/4.),
(**)             ListTools:-MakeUnique(pointList)): 
(**) plots:-display([
(**)     plots:-pointplot(Array(L1),symbol=point,color=red),
(**)     plots:-pointplot(Array(L2),symbol=point),
(**)     plottools:-circle([r/2.,r/2.], r/2.,color=red),
(**)     plottools:-rectangle([0,0],[r,r],color=black, style=line) ],
(**)     scaling=constrained, axes=boxed,
(**)     caption=sprintf("%d point Monte Carlo approximation of Pi = %.4f",
(**)                     nops(pointList), 4. * numC / nops(M)));

Since Monte Carlo methods still produce pretty good answers for all sorts of "bad" PRNGs (for example: the non-prime linear congrence at the bottom of this blog post does just as well), this doesn't seem to be a great test of randomness either. Although, it will almost detect some really egregiously bad generators like the following:

(**) r2a := rand(0..2): r2 := ()-> r2a() mod 2;
(**) L := [seq(r2(), i=1..size*2*1000)]:

Which generates a worse approximation of 2.556 and this plot:

And finally, for something a little bit meta, how about trying to approximate π using the binary digits of π as a source of pseudo-random bits:

(**) approxPi := trunc(evalf[4818+2]((Pi-3)*10^(4818))):
(**) L := convert(approxPi, base, 2): # about 16000 binary digits of Pi

which appears to work about as well as Maple's default Mersenne Twister generated bits:

Here is a worksheet with all the computations from this post: montecarloPi.mw