 Sebastien Binet
Dec 23, 2018 13 min read

# Computing and plotting π with Gonum and a zest of Monte Carlo

Today we will see how we can compute π with a technique called Monte Carlo.

Wikipedia, the ultimate source of truth in the (known) universe has this to say about Monte Carlo:

Monte Carlo methods (or Monte Carlo experiments) are a broad class of computational algorithms that rely on repeated random sampling to obtain numerical results. (…) Monte Carlo methods are mainly used in three distinct problem classes: optimization, numerical integration, and generating draws from a probability distribution.

In other words, the Monte Carlo method is a numerical technique using random numbers.

With Monte Carlo integration, we can estimate the value of an integral:

• take the function value at random points
• the area (or volume) times the average function value estimates the integral.

With Monte Carlo simulation, we can predict an expected measurement:

• an experimental measurement is split into a sequence of random processes
• use random numbers to decide which processes happen
• tabulate the values to estimate the expected probability density function (PDF) for the experiment.

In High Energy Physics (HEP) – but also in many other scientific domains – the Monte Carlo method is used to model a phenomenon, to create a simulation of a given process (and perhaps compare that simulation with measurements of the real world.) In HEP, we have very detailed simulation programs (like Geant4 that models interactions of particles with matter using all our knowledge of particle physics) and fast simulation programs (like Delphes (C++) or fads (in Go)) that very coarsely model physics. But before being able to write a HEP detector simulation, we need to know how to generate random numbers in Go.

## Generating random numbers

The Go standard library provides the building blocks for implementing Monte Carlo techniques, via the math/rand package.

`math/rand` exposes the rand.Rand type, a source of (pseudo) random numbers. With `rand.Rand`, one can:

• generate random numbers following a flat, uniform distribution between `[0, 1)` with `Float32()` or `Float64()`;
• generate random numbers following a standard normal distribution (of mean 0 and standard deviation 1) with `NormFloat64()`;
• and generate random numbers following an exponential distribution with `ExpFloat64`.

If you need other distributions, have a look at Gonum’s gonum/stat/distuv.

`math/rand` exposes convenience functions (`Float32`, `Float64`, `ExpFloat64`, …) that share a global `rand.Rand` value, the “default” source of (pseudo) random numbers. These convenience functions are safe to be used from multiple goroutines concurrently, but this may generate lock contention. It’s probably a good idea in your libraries to not rely on these convenience functions and instead provide a way to use local `rand.Rand` values, especially if you want to be able to change the seed of these `rand.Rand` values.

Let’s see how we can generate random numbers with `"math/rand"`:

 `````` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 `````` ``````package main import ( "fmt" "math/rand" ) func main() { const seed = 12345 src := rand.NewSource(seed) rnd := rand.New(src) const N = 10 for i := 0; i < N; i++ { r := rnd.Float64() // r is in [0.0, 1.0) fmt.Printf("%v\n", r) } }``````

Running this program gives:

 ``` 1 2 3 4 5 6 7 8 9 10 11 ``` ```\$> go run ./mc-0.go 0.8487305991992138 0.6451080292174168 0.7382079884862905 0.31522206779732853 0.057001989921077224 0.9672449323010088 0.6139541710075446 0.01505990819189991 0.13361969083044145 0.5118319569473198```

OK. Does this seem flat to you? Not sure…

Let’s modify our program to better visualize the random data. We’ll use a histogram and the go-hep.org/x/hep/hbook and go-hep.org/x/hep/hplot packages to (respectively) create histograms and display them.

Note: `hplot` is a package built on top of the gonum.org/v1/plot package, but with a few HEP-oriented customization. You can use `gonum.org/v1/plot` directly if you so choose or prefer.

 `````` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 `````` ``````func main() { const seed = 12345 src := rand.NewSource(seed) rnd := rand.New(src) const N = 10000 // create a 1-dim histogram of float64s, with 100 bins, from 0 to 1. huni := hbook.NewH1D(100, 0, 1.0) for i := 0; i < N; i++ { r := rnd.Float64() // r is in [0.0, 1.0) huni.Fill(r, 1) } plot(huni, "uniform.png") }``````

We’ve increased the number of random numbers to generate to get a better idea of how the random number generator behaves, and disabled the printing of the values.

We first create a 1-dimensional histogram `huni` with 100 bins from 0 to 1. Then we fill it with the value `r` and an associated weight (here, the weight is just `1`.)

Finally, we just plot (or rather, save) the histogram into the file `"uniform.png"` with the `plot(...)` function:

 `````` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 `````` ``````func plot(h *hbook.H1D, fname string) { p := hplot.New() // create a new hplot.Plot hh := hplot.NewH1D(h) // create a plotter for the histogram hh.Color = color.NRGBA{0, 0, 255, 255} p.Add(hh, hplot.NewGrid()) const ( width = 10 * vg.Centimeter height = -1 // choose height automatically ) err := p.Save(width, height, fname) if err != nil { log.Fatal(err) } }``````

Running the code creates a `uniform.png` file:

 ```1 ``` `\$> go run ./mc-1.go` Indeed, that looks rather flat.

So far, so good. Let’s add a new distribution: the standard normal distribution.

 `````` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 `````` ``````func main() { const seed = 12345 src := rand.NewSource(seed) rnd := rand.New(src) const N = 10000 huni := hbook.NewH1D(100, 0, 1.0) hgauss := hbook.NewH1D(100, -5, 5) for i := 0; i < N; i++ { r := rnd.Float64() // r is in [0.0, 1.0) huni.Fill(r, 1) g := rnd.NormFloat64() hgauss.Fill(g, 1) } plot(huni, "uniform.png") plot(hgauss, "norm.png") }``````

Running the code creates the following new plot:

 ```1 ``` `\$> go run ./mc-2.go` Note that this has slightly changed the previous `"uniform.png"` plot: we are sharing the source of random numbers between the 2 histograms. The sequence of random numbers is exactly the same than before (modulo the fact that now we generate -at least- twice the number than previously) but they are not associated to the same histograms.

OK, this does generate a gaussian. But what if we want to generate a gaussian with a mean other than `0` and/or a standard deviation other than `1` ?

The math/rand.NormFloat64 documentation kindly tells us how to achieve this:

“To produce a different normal distribution, callers can adjust the output using: `sample = NormFloat64() * desiredStdDev + desiredMean`

Let’s try to generate a gaussian of mean `10` and standard deviation `2`. We’ll have to change a bit the definition of our histogram:

 `````` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 `````` ``````func main() { const seed = 12345 src := rand.NewSource(seed) rnd := rand.New(src) const ( N = 10000 mean = 10.0 stddev = 5.0 ) huni := hbook.NewH1D(100, 0, 1.0) hgauss := hbook.NewH1D(100, -10, 30) for i := 0; i < N; i++ { r := rnd.Float64() // r is in [0.0, 1.0) huni.Fill(r, 1) g := mean + stddev*rnd.NormFloat64() hgauss.Fill(g, 1) } plot(huni, "uniform.png") plot(hgauss, "gauss.png") fmt.Printf("gauss: mean= %v\n", hgauss.XMean()) fmt.Printf("gauss: std-dev= %v +/- %v\n", hgauss.XStdDev(), hgauss.XStdErr()) }``````

Running the program gives:

 ```1 2 3 ``` ```\$> go run mc-3.go gauss: mean= 10.105225624460644 gauss: std-dev= 5.048629091912316 +/- 0.05048629091912316``` Ok, so now we know how to generate random numbers that follow some distribution. How do we evaluate π with that?

## Approximating π using a Monte Carlo technique

Consider a circle inscribed in a unit square:

• the unit square has an area of `d^2 = (2r)^2 = 4r^2`, (where `d` is the diameter of the circle and `r` its radius)
• the unit circle has an area of `π.r^2`.

The ratio of these two areas is thus `area(circle)/area(square) = π/4`.

One can then leverage the Monte Carlo technique to estimate π like so:

• draw a square, inscribe a circle within it,
• uniformly scatter objects of uniform size over the square,
• count the number of objects inside the circle and the total number of objects,
• the ratio of the `inside-count` and the `total-sample-count` is an estimate of the ratio of the two areas, which is `π/4`.

We just have to multiply the result by 4 to estimate π.

We can start with the following code:

 `````` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 `````` ``````package main import ( "flag" "fmt" ) func main() { n := flag.Int("n", 1e7, "MC sample size") flag.Parse() fmt.Printf("pi(%d) = %1.16f\n", *n, pi(*n)) } func pi(n int) float64 { // ??? }``````

We just have to fill in the blanks :).

Following the algorithm laid out above, we know we have to draw `n` objects randomly over the unit square and count the number of objects (or darts if you are into this kind of game) that fall inside the circle. An object can be identified by its `(x,y)` coordinates: we have to draw 2 random values `x` and `y` between `0` and `1` (the dimensions of the top-right quarter-square of the unit square.) This can be translated into the following code:

 ``````1 2 3 4 `````` ``````import "math/rand" x := rand.Float64() y := rand.Float64()``````

Then, deciding whether this `(x,y)` dart is inside the (quarter) circle is just a matter of applying Pythagoras’ theorem: `x*x + y*y < 1`.

i.e.:

 `````` 1 2 3 4 5 6 7 8 9 10 11 `````` ``````func pi(n int) float64 { inside := 0 for i := 0; i < n; i++ { x := rand.Float64() y := rand.Float64() if (x*x + y*y) < 1 { inside++ } } return 4 * float64(inside) / float64(n) }``````

Let’s Go!

 ``` 1 2 3 4 5 6 7 8 9 10 ``` ```\$> for i in `seq 1 9`; do go run ./mc-pi.go -n=`echo '10^'\$i | bc`; done pi(10) = 3.6000000000000001 pi(100) = 3.3599999999999999 pi(1000) = 3.1680000000000001 pi(10000) = 3.1600000000000001 pi(100000) = 3.1518000000000002 pi(1000000) = 3.1405520000000000 pi(10000000) = 3.1414072000000002 pi(100000000) = 3.1415181200000002 pi(1000000000) = 3.1415852000000002```

Ok… Sadly, this is perhaps not a very quickly converging method…

### Graphics

Just for fun, let’s add a little GUI part to visualize where the darts land.

Our GUI will be a web server with two end points:

• `/` will plot a quarter circle inscribed inside the top-right quarter unit square,
• `/data` will send plots over a WebSocket.

 `````` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 `````` ``````func main() { log.SetPrefix("monte-carlo: ") log.SetFlags(0) n := flag.Int("n", 1e7, "number of samples") flag.Parse() srv := newServer() go srv.pi(*n) http.HandleFunc("/", plotHandle) http.Handle("/data", websocket.Handler(srv.dataHandler)) log.Printf("listening on :8080...") err := http.ListenAndServe(":8080", nil) if err != nil { log.Fatal(err) } }``````

We will start with creating a web server:

 ``````1 2 3 4 5 6 7 8 `````` ``````type server struct { in plotter.XYs // points inside the circle out plotter.XYs // points outside the circle n int // number of samples datac chan float64 // channel of (x,y) points randomly drawn plots chan wplot // channel of base64-encoded PNG plots }``````

The `in` and `out` fields are slices of `(x,y)` points that implement the `plotter.XYer` interface:

 ```1 2 3 4 5 6 7 ``` ```\$> go doc gonum.org/v1/plot/plotter.XYs type XYs []struct{ X, Y float64 } XYs implements the XYer interface. func (xys XYs) Len() int func (xys XYs) XY(i int) (float64, float64)```
 ```1 2 3 4 5 6 7 8 9 ``` ```\$> go doc gonum.org/v1/plot/plotter.XYer type XYer interface { // Len returns the number of x, y pairs. Len() int // XY returns an x, y pair. XY(int) (x, y float64) } XYer wraps the Len and XY methods.```

The `plotter.XYer` interface is used by gonum/plot to plot `(x,y)` points.

Values of type `server` will be created via the `newServer` function:

 `````` 1 2 3 4 5 6 7 8 9 10 11 12 `````` ``````func newServer() *server { srv := &server{ in: make(plotter.XYs, 0, 1024), out: make(plotter.XYs, 0, 1024), datac: make(chan float64), plots: make(chan wplot), } go srv.run() return srv }``````

The `run()` method is a simple `for`-loop that listens on the `srv.datac` channel and sends the resulting plot on the `srv.plots` channel. The `srv.datac` channel is filled in the `pi()` method:

 ``````1 2 3 4 5 6 7 `````` ``````func (srv *server) pi(samples int) { for i := 0; i < samples; i++ { x := rand.Float64() y := rand.Float64() srv.datac <- float64{x, y} } }``````
 `````` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 `````` ``````type Point struct { X, Y float64 } func (srv *server) run() { for v := range srv.datac { srv.n++ x := v y := v d2 := x*x + y*y pt := Point{x, y} switch { case d2 < 1: srv.in = append(srv.in, pt) default: srv.out = append(srv.out, pt) } srv.plots <- plot(srv.n, srv.in, srv.out) } }``````

Next is the `plot` function:

 `````` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 `````` ``````func plot(n int, in, out plotter.XYs) wplot { radius := vg.Points(0.1) p := hplot.New() p.X.Label.Text = "x" p.X.Min = 0 p.X.Max = 1 p.Y.Label.Text = "y" p.Y.Min = 0 p.Y.Max = 1 pi := 4 * float64(len(in)) / float64(n) p.Title.Text = fmt.Sprintf("n = %d\nπ = %v", n, pi) sin, err := hplot.NewScatter(in) if err != nil { log.Fatal(err) } sin.Color = color.RGBA{255, 0, 0, 255} // red sin.Radius = radius sout, err := hplot.NewScatter(out) if err != nil { log.Fatal(err) } sout.Color = color.RGBA{0, 0, 255, 255} // blue sout.Radius = radius p.Add(sin, sout, hplot.NewGrid()) return wplot{Plot: renderImg(p)} }``````

We create a new `hplot.Plot` – a thin wrapper around the `plot.Plot` type from gonum/plot – with the correct labels. We plot the points inside the circle in red and the ones outside in blue.

The `wplot` type is just a shim to hold the resulting `base64` encoded string of the PNG plot, created with `renderImg`:

 ``````1 2 3 `````` ``````type wplot struct { Plot string `json:"plot"` }``````
 `````` 1 2 3 4 5 6 7 8 9 10 11 `````` ``````func renderImg(p *hplot.Plot) string { size := 20 * vg.Centimeter canvas := vgimg.PngCanvas{vgimg.New(size, size)} p.Draw(draw.New(canvas)) out := new(bytes.Buffer) _, err := canvas.WriteTo(out) if err != nil { log.Fatal(err) } return base64.StdEncoding.EncodeToString(out.Bytes()) }``````

At the other end of the `srv.plots` channel is the `dataHandler` method that pulls plots out to send them to the web client, over a WebSocket:

 ``````1 2 3 4 5 6 7 8 `````` ``````func (srv *server) dataHandler(ws *websocket.Conn) { for data := range srv.plots { err := websocket.JSON.Send(ws, data) if err != nil { log.Printf("error sending data: %v\n", err) } } }``````

and, finally, the `/` end point and its home page:

 ``````1 2 3 `````` ``````func plotHandle(w http.ResponseWriter, r *http.Request) { fmt.Fprintf(w, page) }``````
 `````` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 `````` ``````const page = ` Monte Carlo

```````

Running the Go program and navigating to `localhost:8080`, you should see:

 ```1 2 ``` ```\$> go run ./mc-pi-plot.go monte-carlo: listening on :8080...``` ## Conclusions

I hope this quick foray into a technique that is at the heart of many physics simulations was fun. The Monte Carlo technique isn’t always the fastest technique to perform simulations (this obviously depends on the “shape” of the function we want to model) but it is deceptively simple, and can be visually “fun” – for some definition of fun.

Note that in this little example, we had fun creating a PNG image, encoding it to JSON+`base64` and sending it over a WebSocket. This was just to exercize a bunch of packages with Go. Gonum/plot has support for a few backends in addition to PNG: JPEG, TIFF, EPS, PDF, SVG, TeX.

The complete code for `mc-pi-plot.go` is here.