library(dplyr)
library(ggplot2)
library(knitr)
opts_chunk$set(message = FALSE, eval = TRUE)
Simulate Buffon’s Needle
n <- 20
d <- runif(n = n, min = 0, max = 1/2)
theta <- runif(n = n, min = 0, max = pi)
(df <- data_frame(d, theta))
## # A tibble: 20 x 2
## d theta
## <dbl> <dbl>
## 1 0.421537737 0.08101089
## 2 0.080696171 1.10375158
## 3 0.008729169 1.19186020
## 4 0.449464490 0.83817382
## 5 0.328243365 0.88036281
## 6 0.499262640 0.53429573
## 7 0.372272875 0.07884363
## 8 0.498089879 0.73014834
## 9 0.178281839 2.28789951
## 10 0.157958551 1.30119270
## 11 0.466539003 1.23536838
## 12 0.402765643 1.05361257
## 13 0.474079798 2.25967411
## 14 0.366149133 2.10208444
## 15 0.393188542 1.05431343
## 16 0.439498770 3.13983823
## 17 0.244888635 1.29131678
## 18 0.197206581 0.55139906
## 19 0.080525428 1.37555721
## 20 0.327939279 0.58738903
sin(theta)/2 > d
## [1] FALSE TRUE TRUE FALSE TRUE FALSE FALSE FALSE TRUE TRUE TRUE
## [12] TRUE FALSE TRUE TRUE FALSE TRUE TRUE TRUE FALSE
mean(sin(theta)/2 > d)
## [1] 0.6
2/pi
## [1] 0.6366198
Visualize
(p <- ggplot(data = df, aes(x = theta, y = d)) +
scale_x_continuous(expand = c(0, 0), limits = c(0, pi)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 1/2)) +
stat_function(fun = function(x) {sin(x)/2}) +
theme_bw())
df <- mutate(df, cross = sin(theta)/2 > d)
p + geom_point(data = df, aes(color = cross))
Estimating Pi
2 * n / sum(df$cross)
## [1] 3.333333