Finding #pi with #montecarlo method and #Clojure — #math #justmath

I was reading a post from Toward Data Science blog this morning on mathematical programming to build up skills in data science by . While the article was based around Python it didn’t use any of the popular frameworks like NumPy or SciPy.

Now with a bit of a lull I wanted to keep my brain ticking nicely so the thought of using math within Clojure appeals nicely to me. And I’m not saying one is better than the other. The best language to for data science is the one you know. The main key of data science is having a good grounding of the math behind, not the frameworks to make it easier.

Calculating Pi By Simulating Random Dart Board Throws

The Monte Carlo method is the concept of emulating a random process. When the process is repeated a large number of times will give rise to the approximation of some mathematical quantity of interest.

If you imagine a square dart board…..

Now imagine a square dart board with a circle inside the square, the edges of circle touch the square…..

If you throw enough darts at the board some will land within the circle and some outside of it. As the original article graphically put it:

These are random throws, you might throw 10 times, you might throw 1 million times. At the end of the dart throws you count the number of darts within the circle, divide that by the number of throws (10, 1m etc) and then multiply it by 4.

As the original article states: the probability of a dart falling inside the circle is just the ratio of the area of the circle to that of the area of the square board.

The more throws we do the better chance we get of finding a number near Pi. The law of large numbers at work.

Throwing a Dart at The Board

I’m going to create a function that simulates a single dart throw. I want to break down my Clojure code into as many simple functions as possible. This makes testing and bug finding far easier in my opinion.

(defn throw-dart [] {:x (calc-position 0) :y (calc-position 0)})

What I’m creating an x,y coordinate with a 0,0 centre point then passing the coord for the x and the y through another function to calculate the position ( calc-position).

(def side-of-square 2) (defn calc-position [v] (* (/ (+ v side-of-square) 2) (+ (- 1) (* 2 (Math/random)))))

The calc-position function takes the value of either x and y and applies the calculation, this is somewhere -side-of-square/2 and +side-of-square/2 around the centre point.

Running this function in a REPL we can see the x or y positions.

mathematical.programming.examples.montecarlo> (calc-position 0) 0.4298901518005238

Is The Dart Within The Circle?

Now I have a x,y position as a map {:x some random throw value :y some random throw value} I want to confirm that the throw is within the circle.

Using the side-of-square value again (hence it's a def ) I can figure out if the dart hits within. I'll pass the map with x,y coords in and take the square root of the added squared coordinates.

(defn is-within-circle [m] (let [distance-from-center (Math/sqrt (+ (Math/pow (:x m) 2) (Math/pow (:y m) 2)))] (< distance-from-center (/ side-of-square 2))))

This function will return true or false. If I check this in the REPL it looks like this:

mathematical.programming.examples.montecarlo> (throw-dart) {:x 0.22535085231582297, :y 0.04203583357796781} mathematical.programming.examples.montecarlo> (is-within-circle *1) true

Now Throws Lots of Darts

So far there are functions to simulate a dart throw and confirm it’s within the circle. Now I need to repeat this process as many times as required.

I’m creating two functions, compute-pi-throwing-dart to run a desired number of throws and throw-range to do the actual working to find the number of true hits in the circle.

(defn throw-range [throws] (filter (fn [t] (is-within-circle (throw-dart))) (range 0 throws))) (defn compute-pi-throwing-dart [throws] (double (* 4 (/ (count (throw-range throws)) throws))))

The throw-range function executes the throw-dart function and is-within-circle evaluates the map to see if the value is either true or false. The filter functions will return a list of true values. So, for example, if out of ten throws the first, third and fifth are within the circle I'll get (1,3,5) as the result from the function.

Calling the function compute-pi-throwing-dart sets all this into motion. Like I said at the start, taking the number of darts in the circle and dividing that by the number of throws taken, multiplying that by four should give a number close to Pi.

The more throws you do, the closer it should get.

mathematical.programming.examples.montecarlo> (compute_pi_throwing_dart 10) 3.2 mathematical.programming.examples.montecarlo> (compute_pi_throwing_dart 10) 3.2 mathematical.programming.examples.montecarlo> (compute_pi_throwing_dart 10) 3.6 mathematical.programming.examples.montecarlo> (compute_pi_throwing_dart 10) 2.4 mathematical.programming.examples.montecarlo> (compute_pi_throwing_dart 10) 4.0 mathematical.programming.examples.montecarlo> (compute_pi_throwing_dart 10) 2.8 mathematical.programming.examples.montecarlo> (compute_pi_throwing_dart 100) 2.92 mathematical.programming.examples.montecarlo> (compute_pi_throwing_dart 1000) 3.136 mathematical.programming.examples.montecarlo> (compute_pi_throwing_dart 10000) 3.138 mathematical.programming.examples.montecarlo> (compute_pi_throwing_dart 100000) 3.15456 mathematical.programming.examples.montecarlo> (compute_pi_throwing_dart 1000000) 3.13834 mathematical.programming.examples.montecarlo> (compute_pi_throwing_dart 10000000) 3.1419096

Let’s Build a Simulation

Via the REPL there is proof of an emergent behaviour, the value of Pi comes from the large number of throws we did at the dart board.

The last thing I’ll do is build a function to run the simulation.

(defn run-simulation [iter] (map (fn [i] (let [throws (long (Math/pow 10 i))] (compute-pi-throwing-dart throws))) (range 0 iter)))

If I run 4 simulations I’ll get 1, 10, 100 and 1000 throws computed, these are then returned as a list. If I run 9 simulations (which can take some time depending on the machine you’re using) in the REPL I get the following:

mathematical.programming.examples.montecarlo> (run-simulation 9) (0.0 3.6 3.28 3.128 3.1176 3.1428 3.142932 3.1425368 3.14173752)

That’s a nice approximation, Pi is 3.14159265 so to get a Monte Carlo method to compute Pi by random evaluations is good.

Originally published at http://dataissexy.wordpress.com on April 15, 2019.