Suppose that you drop a short needle on ruled paper, what is then the probability that the needle comes to lie in a position where it crosses one of the lines? - Georges Buffon, 1777
When the needle is the right length, the answer is pi!
Here is a quick example in R of approximating the value of Pi by digitally tossing a box of matches onto a sheet of lined paper.
library(retistruct)Make a needle
match_stick_length <- 20
number_of_matches_in_a_box <- 10000
distance_between_surface_lines <- 40
point_1 <- c(-(match_stick_length/2),(match_stick_length/2))
point_2 <- c(0,0)
single_match <-cbind(point_1, point_2)
rownames(single_match)<- c("point_1","point_2")
colnames(single_match)<- c("x","y")
single_match## x y
## point_1 -10 0
## point_2 10 0## Warning: package 'rgl' was built under R version 3.4.4## Warning: package 'shiny' was built under R version 3.4.4
Throw the match
First rotate around the origin
angle <- runif(1, min=1, max=360)*pi/180
xy <- as.matrix(single_match)
# Rotation
cos.angle <- cos(angle)
sin.angle <- sin(angle)
after_rotation <- xy %*% t(matrix(c(cos.angle, sin.angle, -sin.angle,
cos.angle), 2, 2))
after_rotation## [,1] [,2]
## point_1 0.1190636 -9.999291
## point_2 -0.1190636 9.999291
Then translate that line onto the table
# Translation
translation.x <- runif(1, min=10, max=90)
translation.y <- runif(1, min=10, max=90)
translation <- matrix(c(translation.x,translation.x,translation.y,translation.y), 2, 2)
new_coord <- after_rotation + translation
colnames(new_coord) <- c("x","y")
rownames(new_coord) <- c("Point_1","Point_2")new_coord <- after_rotation + matrix(c(50,50,50,50),2,2)
Repeat for an entire box of matches
Throwing_matches <- function(){
single_match <-cbind(c(-5,5),c(0,0))
angle <- runif(1, min=1, max=360)*pi/180
xy <- as.matrix(single_match)
# Rotation
cos.angle <- cos(angle)
sin.angle <- sin(angle)
after_rotation <- xy %*% matrix(c(cos.angle, sin.angle, -sin.angle,
cos.angle), 2, 2, byrow=TRUE )
# Translation
translation.x <- runif(1, min=10, max=90)
translation.y <- runif(1, min=10, max=90)
translation <- matrix(c(translation.x,translation.x,translation.y,translation.y), 2, 2)
new_coord <- after_rotation + translation
colnames(new_coord) <- c("x","y")
rownames(new_coord) <- c("Point_1","Point_2")
return(new_coord)
}
Throwing_matches()## x y
## Point_1 50.61685 36.22574
## Point_2 41.00933 38.99984for(i in 1:1000){
Throwing_matches()
}
surface_line <- function(line_height = 0){matrix(c(0,100,line_height,line_height), 2, 2, byrow=FALSE)}## [,1] [,2]
## [1,] 0 20
## [2,] 100 20
one_match <- Throwing_matches()
one_surface_line <- surface_line(30)
line.line.intersection(one_match[1,], one_match[2,], one_surface_line[1,], one_surface_line[2,], interior.only = TRUE)## [1] NA NA
Make a list of surface lines and a list of matches
list_of_matches <- array(NA,dim = c(2,2,number_of_matches_in_a_box))
list_of_surface_lines <- array(NA,dim = c(2,2,6))
for(i in 1:number_of_matches_in_a_box){
list_of_matches[,,i] <- Throwing_matches()
}
line_placement <- seq(0,100, by=20)
for(i in 1:6){
list_of_surface_lines[,,i] <- surface_line(line_placement[i])
}
crosses <- matrix(NA, 1, 10, byrow=TRUE)
for(i in 1:dim(list_of_matches)[3]){
for(j in 1:dim(list_of_surface_lines)[3]){
crosses <- rbind(crosses, c(
line.line.intersection(list_of_matches[1,,i], list_of_matches[2,,i], list_of_surface_lines[1,,j],
list_of_surface_lines[2,,j], interior.only = TRUE),
list_of_matches[1,,i], list_of_matches[2,,i], list_of_surface_lines[1,,j], list_of_surface_lines[2,,j])
)
}
}
crosses <- crosses[which(is.na(crosses[,1]) == FALSE),]
colnames(crosses) <- c("x.cross", "y.cross", "x.p1", "y.p1", "x.p2",
"y.p2", "x.surface.p1", "y.surface.p1", "x.surface.p2",
"y.surface.p2")
head(crosses)## x.cross y.cross x.p1 y.p1 x.p2 y.p2 x.surface.p1
## [1,] 12.31511 20 14.36624 22.62863 8.214425 14.74477 0
## [2,] 45.79014 80 44.77192 82.43566 48.628934 73.20942 0
## [3,] 70.15395 20 74.24544 28.20362 69.782313 19.25485 0
## [4,] 49.66294 60 49.75337 59.89791 43.122492 67.38333 0
## [5,] 59.91925 20 60.98953 22.25686 56.704609 13.22141 0
## [6,] 74.89057 60 72.66679 55.32310 76.960902 64.35419 0
## y.surface.p1 x.surface.p2 y.surface.p2
## [1,] 20 100 20
## [2,] 80 100 80
## [3,] 20 100 20
## [4,] 60 100 60
## [5,] 20 100 20
## [6,] 60 100 60
pi_estimate <- (2* match_stick_length * number_of_matches_in_a_box)/(distance_between_surface_lines * length(crosses[,1]))## [1] 3.194888