Introduction

I knew of Monte Carlo simulations in my work with planning cellular networks using radio planning software. These software tools are specialised systems used to aid cellular network design, site selection and capacity modelling. The latter is often implemented using Monte Carlo simulations to model the behaviour of users, their geographic locations within the radio cell and their services. Simplistically speaking, when using these tools you generally tick a few boxes, enter some values and hit ‘run’. Whilst I have used a number of radio planning systems, my knowledge of how these statistical models worked below the surface was sketchy. I could tell someone what a Monte Carlo simulation was used for but I had no idea how it worked, so I thought it might be interesting to try some examples using R.

Definition

Monte Carlo methods cover a broad range of applications and refer to the use of algorithms employing random samples to derive a numerical solution for deterministic problems (i.e. those in which randomness plays no role). These methods can be used to derive the outcome of any problem that has uncertainty in terms of its input variables. By the law of large numbers the average of the results of repeated simulations should converge on the expected value as the number of iterations is increased.

I read a number of definitions online and in several textbooks and frankly I wasn’t much the wiser on the topic. It wasn’t until I played around with a few examples that it made a little more sense. The example provided on the main wikipedia entry seemed like a fairly interesting application, so I tried my hand at reproducing it using R.

Calculating the Value of Pi

In this example, a Monte Carlo simulation is used to calculate the value of pi, using simple geometry and randomly generated points. Consider a unit square (i.e. a square with an area equal to one) containing a quadrant of a circle with radius equal to 1 and an area equal to $\frac{\pi.r^2}{4} = \frac{\pi}{4}$. It follows that the ratio of their areas is thus $\frac{\pi}{4}$.

There are several ways of plotting a circle using R and I chose to do it using some familiar equations and plotly. First I created a sequence for the angles in the quadrant, and then calculated the values of $x$ and $y$.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

# Define angle from 0 to 90 degrees 
theta <- seq(0, 90, 0.1)

# Convert to radians
theta <- theta * pi/180

# Circle parameters 
r <- 1
x <- r * cos(theta)
y <- r * sin(theta)

# Create data frame
circ <- data.frame(x, y)

Once I had created the data frame, the plot was simple enough. The circle is drawn using a scatter plot and the square is added in using the shapes() function with appropriate boundary values.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

# Plot
plot_ly(circ,
        x = ~x,
        y = ~y,
        type= 'scatter',
        mode = 'markers',
        width = 600, 
        height = 600) %>%
layout(showlegend = FALSE, 
       shapes = list(type = "rect",
                     line = list(color = "blue"),
                     opacity = 0.3,
                     x0 = 0, x1 = 1, xref = "x",
                     y0 = 0, y1 = 1, yref = "y")) %>%
config(displayModeBar = F) # Hide the mode bar

If points are uniformly plotted within the unit square, then the proportion of the points inside the quadrant compared to the total number of points allows one to approximate the value of $\pi$, by observing that:

$$\frac{\Sigma points_{inside}}{\Sigma points_{total}} \approx \frac{\pi}{4}$$
Distributing the points randomly within the unit square is where the Monte Carlo method comes into play. Points are allocated an $(x,y)$ coordinate with $0 < x \leq 1$ and $0 < y \leq 1$. In R, this can be done using the runif() function. Points inside the quadrant are coloured blue and those outside are coloured red. The approximate value of $\pi$ can then be determined for the number of samples $n$.

 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
36
37
38
39
40
41
42
43
44

# Number of samples
n = 1000

# x coordinates
ptx <- runif(n)
# y coordinates
pty <- runif(n)

# Make a data frame
dfpt <- data.frame(ptx, pty)
  
# Calculate the radius (Pythagoras)
dfpt$ptr <- sqrt(dfpt$ptx^2 + dfpt$pty^2)
  
# Add in color of point
# If its inside then blue, otherwise red
dfpt$col <- ifelse(dfpt$ptr < 1, "blue", "red")

# Calculate PI
pt_in <- length(which(dfpt$col == "blue"))
pt_out <- length(which(dfpt$col == "red"))
dfpt$pi <- pt_in/(pt_out+pt_in)*4

# Plotly
plot_ly(dfpt, # circle
        x = ~x,
        y = ~y,
        type= 'scatter',
        mode = 'markers',
        width = 600, 
        height = 600) %>%

layout(showlegend = FALSE,
       shapes = list(type = "rect", # square
                     line = list(color = "blue"),
                     opacity = 0.3,
                     x0 = 0, x1 = 1, xref = "x",
                     y0 = 0, y1 = 1, yref = "y")) %>%
  
add_markers(x = ~dfpt$ptx, # data points
            y = ~dfpt$pty,
            marker = list(color = ~dfpt$col))  %>%

config(displayModeBar = F)

Even with only 1000 points, the approximation of $\pi$ is correct to two significant digits. Increasing the number of samples $n$ generally improves the accuracy of the estimation, but in some cases the individual results may be worse. The real value in this method is to run the simulation multiple times and then average the results.

 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

# Number of samples
n = 1000

# Number of simulations
s = 100

# Initialise results object
pi_result <- c()

# Loop for the number of simulations
for (i in 1:s){
  
  # Assign random coordinates
  # x coordinates
  ptx <- runif(n)
  # y coordinates
  pty <- runif(n)
  
  # Coerce
  dfpt <- data.frame(ptx, pty)
  
  # Add in radius
  dfpt$ptr <- sqrt(dfpt$ptx^2 + dfpt$pty^2)
  
  # Add in color of point
  dfpt$col <- ifelse(dfpt$ptr < 1, "blue", "red")

  # Calculate PI
  pt_in <- length(which(dfpt$col == "blue")) # Points inside
  pt_out <- length(which(dfpt$col == "red")) # Points outside
  pi_calc <- pt_in/(pt_out+pt_in)*4
  pi_tmp <- data.frame(i, pi_calc)
  pi_result <- rbind(pi_result, pi_tmp)
}

Running the simulation 100 times and plotting the result for each simulation shows how the estimates can vary considerably. I tabulated the results for the first 10 simulations to demonstrate this effect.

 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

# Mean result
pi_result$pi_mean <- mean(pi_result$pi_calc)
# Pi value to 6 sig. digits
pi_result$pi_actual <- 3.14159

# Legend
l <- list(bordercolor = "#D3D3D3",
      borderwidth = 2,
      xanchor = 'center',
      orientation = 'h',
      y = 1.1,
      x = 0.5) 

# Plot
plot_ly(pi_result,
        x = ~i,
        y = ~pi_calc,
        name = 'Estimate',
        type = 'scatter',
        mode = 'markers') %>%
add_lines(pi_result, y = ~pi_actual, name = 'Pi') %>% 
add_lines(pi_result, y = ~pi_mean, name = 'Mean Estimate') %>%
config(displayModeBar = F) %>%
layout(yaxis = list(title = "Pi"),
       xaxis = list(title = "Simulation Number"),
       legend = l)

In the preceding plot, the mean estimate and the value of $\pi$ overlap, so plotting the rolling mean of the estimations is a bit more helpful. It shows how the approximation converges towards the value of $\pi$ as the simulation count increases.

Finally, I set $n$ = 10,000 and simulations $s$ = 10,000 to see how this improved the result. The final approximation was 3.141523, which is accurate to 5 significant digits.

Conclusion

This post has provided a simple example of how to perform and visualise a Monte Carlo simulation using R. Monte Carlo simulations are used in many domains in engineering and science and I hope to revisit this with some other applications.

Other Reading

Wikipedia: https://en.wikipedia.org/wiki/Monte_Carlo_method

Multiple Examples: https://www.countbayesie.com/blog/2015/3/3/6-amazing-trick-with-monte-carlo-simulations

Manufacturing Example: https://www.programmingr.com/monte-carlo-simulation-in-r/