I like to think of Monte Carlo as a counting method. If a condition is satisfied we make a note (e.g., 1), and if the condition is not satisfied we make a different note (e.g., 0). We then iterate and evaluate the pattern of 1’s and 0’s to learn about our process. Art can be described in a similar way: if a condition is satisfied we use a color, and if a condition is not satisfied we use a different color. After many iterations, we have an image.

Here is a simulation that “draws” a process, inspired by Caleb Madrigal (link here).

# The Data Generating Process

```
f <- function(x){
2*sin(4*x) + 2*sin(5*x) + 12
}
```

# Some Initial Values

`x <- seq(0, 10, length.out = 1000)`

# Using the DGP to generate values of Y

```
y <- f(x)
plot(x, y)
```

# This is the process we want to “draw”

# Now for the Monte Carlo

We are going to evaluate 10,000 points within our process space (10 x 16).

```
num_points <- 10000
rect_width <- 10
rect_height <- 16
points <- matrix(, ncol = 2, nrow = num_points)
```

Column 1 of our points matrix represents the width of our process space while column 2 represents its height. First we fill the matrix with random values within our process space:

```
for(i in 1:num_points){
points[i,1] = runif(1, 0, rect_width)
points[i,2] = runif(1, 0, rect_height)
}
```

Now we iterate across all of those points and evaluate them with respect to our process. Think of the “width” as *X* values and the “height” as *Y* values. Given a value of *X*, is our random value of *Y* less than it would be if we created a *Y* value by using our function (f(x))? If so, mark it in the “points_under” vector. If not, mark it in the “points_over” vector.

```
points_under = matrix(, ncol = 2, nrow = num_points)
points_above = matrix(, ncol = 2, nrow = num_points)
for(i in 1:num_points){
if(points[i,2] < f(points[i,1])){
points_under[i,1] <- points[i,1]
points_under[i,2] <- points[i,2]
}
else{
points_above[i,1] <- points[i,1]
points_above[i,2] <- points[i,2]
}
}
```

Put the results into new vectors without NA’s. Some NA’s come up because our data generating process is crazy.

```
points_under_x <- points_under[!is.na(points_under[,1]),1]
points_under_y <- points_under[!is.na(points_under[,2]),2]
points_over_x <- points_above[!is.na(points_above[,1]),1]
points_over_y <- points_above[!is.na(points_above[,2]),2]
```

Now we have an image…

`plot(points_under_y ~ points_under_x, pch = 20, cex = 0.3)`

Bo\(^2\)m =)