http://web.mit.edu/tee/www/bertrand/problem.html

bertrands-paradox

suppressPackageStartupMessages({                                           
  library(tidyverse)                                                         
  library(sf)                                                                
})                                                                         

# c.f. http://web.mit.edu/tee/www/bertrand/problem.html                    

chord_length <- function(theta) {                                          
  2 * sin(theta / 2)                                                         
}                                                                          

point <- function(x, y) {                                                  
  list(x = x, y = y)                                                         
}                                                                          

angle_to_point <- function(angle) {                                        
  point(x = cos(angle),                                                      
        y = sin(angle))                                                            
}                                                                          

euclidean_distance <- function(a, b) {                                     
  sqrt((a$y - b$y)^2 + (a$x - b$x)^2)                                        
}                                                                          

# Picking random angles                                                    
runif(1e4L, min = 0, max = 2 * pi) %>%                                     
  map_dbl(chord_length) %>%                                                  
  map_lgl(~ . > sqrt(3)) %>%                                                 
  mean()                                                                     
#> [1] 0.3354

# Picking any two random points on the circle                              
tibble(a1 = runif(1e4L, min = 0, max = 2 * pi),                            
       a2 = runif(1e4L, min = 0, max = 2 * pi)) %>%                               
  mutate(p1 = map(a1, angle_to_point),                                       
         p2 = map(a2, angle_to_point),                                              
         d  = map2_dbl(p1, p2, euclidean_distance)) %>%                             
  pull(d) %>%                                                                
  map_lgl(~ . > sqrt(3)) %>%                                                 
  mean()                                                                     
#> [1] 0.3253

# Pick a random point within the circle, and make a chord with that as it's
# midpoint                                                                 
distance_from_midpoint <- function(x, y) {                                 
  # Create the unit circle                                                   
  theta <- seq(0, 2 * pi, 0.001 * pi)                                        
  circle <-                                                                  
    cbind(x = cos(theta),                                                      
          y = sin(theta)) %>%                                                        
    st_linestring()                                                            
  
  # Create the chord                                                         
  chord_gradient <- -1 * (y / x)                                             
  chord_intercept <- y - chord_gradient*x                                    
  xs <- seq(-5, 5, 0.001)                                                    
  chord <-                                                                   
    cbind(x = xs,                                                              
          y = chord_gradient * xs + chord_intercept) %>%                             
    st_linestring()                                                            
  
  points <-                                                                  
    st_intersection(circle, chord) %>%                                         
    as.matrix() %>%                                                            
    as_tibble() %>%                                                            
    `colnames<-`(c("x", "y")) %>%                                              
    pmap(point)                                                                
  
  euclidean_distance(points[[1]], points[[2]])                               
}                                                                          

tibble(r     = runif(1e4L),                                                
       theta = runif(1e4L, min = 0, max = 2 * pi)) %>%                            
  mutate(x = sqrt(r) * cos(theta),                                           
         y = sqrt(r) * sin(theta),                                                  
         d = map2_dbl(x, y, distance_from_midpoint)) %>%                            
  pull(d) %>%                                                                
  map_lgl(~ . > sqrt(3)) %>%                                                 
  mean()                                                                     
#> [1] 0.6087