5.10 Random Data Generation

PyXPlot has a range of mathematical functions which draw random samples from a variety of probability distributions. These are:

These functions all rely upon a common underlying random number generator, whose seed may be set using the set seed command, which should be followed by any integer.

Using random numbers to estimate the value of $\pi $.

PyXPlot’s functions for generating random numbers are most commonly used for adding noise to artificially-generated data. In this example, however, we use them to implement a rather inefficient algorithm for estimating the value of the mathematical constant $\pi $. The algorithm works by spreading randomly-placed samples in the square $\left\{  -1<x<1;\;  -1<y<1\right\} $. The number of these which lie within the circle of unit radius about the origin are then counted. Since the square has an area of $4\, \mathrm{unit}^2$ and the circle an area of $\pi \, \mathrm{unit}^2$, the fraction of the points which lie within the unit circle equals the ratio of these two areas: $\pi /4$.

The following script performs this calculation using $N=5000$ randomly placed samples. Firstly, the positions of the random samples are generated using the random() function, and written to a file called random.dat using the tabulate command. Then, the foreach datum command – which will be introduced in Section 6.5 – is used to loop over these, counting how many lie within the unit circle.

Nsamples = 5000

set samp Nsamples
set output "random.dat"
tabulate 1-2*random():1-2*random() using 0:2:3

n=0
foreach datum i,j in ’random.dat’ using 2:3 {
n = n + (hypot(i,j)<1)
}

print "pi=%s"%(n / Nsamples * 4)

On the author’s machine, this script returns a value of $3.1352$ when executed using the random samples which are returned immediately after starting PyXPlot. This method of estimating $\pi $ is well modelled as a Poisson process, and the uncertainty in this result can be estimated from the Poisson distribution to be $1/\sqrt {N}$. In this case, the uncertainty is $0.01$, in close agreement with the deviation of the returned value of $3.1352$ from more accurate measures of $\pi $.

With a little modification, this script can be adapted to produce a diagram of the datapoints used in its calculation. Below is a modified version of the second half of the script, which loops over the datapoints stored in the datafile random.dat. It uses PyXPlot’s vector graphics commands, which will be introduced in Chapter 3, to produce such a diagram:

set multiplot ; set nodisplay
# Draw a unit circle and a unit square
box from -width/2,-width/2 to width/2,width/2
circle at 0,0 radius width/2 with lt 2

# Now plot the positions of the random data points
n=0
foreach datum i,j in ’random.dat’ using 2:3
{
point at width/2*i , width/2*j with ps 0.1
n = n + (hypot(i,j)<1)
}
set display ; refresh

print "pi=%s"%(n / Nsamples * 4)

The graphical output from this script is shown below. The number of datapoints has been reduced to Nsamples$=500$ for clarity:

\includegraphics{examples/eps/ex_pi_estimation}