Intro

Bertrand’s Paradox 是在機率論裡頭一個著名的例子，旨在表達若我們沒有清楚定義是怎麼樣的「隨機」時，很可能會對機率的詮釋出現 ambiguity。

這個悖論的描述請參考維基百科 ，大意是：「有一圓，內接一正三角形。問隨機地取出一條弦，此弦長大於內接正三角形邊長的機率為何？」

維基中也描述了三種 approach來取弦，分別有：

1. random endpoints：在圓上任取兩點，構成一弦
2. random midpoint：在圓內任取一點，找出以此點為中點的弦。
3. random radial point：在圓上任取一點，與圓心連線，在此線段上任取一點作垂線，構成一弦。

看到了這部影片後，我就想試著用 R 來模擬看看，並繪製 gif 動畫

Animation

Approach 1

第一種方法得到的 empirical probability 約是 1/3

View post on imgur.com

Approach 2

第二種方法得到的 empirical probability 約是 1/4。可以發現靠近圓心的區域比較稀疏。

View post on imgur.com

Approach 3

第三種方法得到的 empirical probability 約是 1/2

View post on imgur.com

小結

這邊想引用 3Blue1Brown 頻道作者的話 ”You have to be careful when the probabilities with an infinite space are involved”，實際上我們可以從模擬的過程發現，針對三種作法，我們「取隨機數」的方式是不一樣的。

針對第一種方法，我們是「在圓上取隨機數」，可以想像成我們用某種方式把圓剪開來，攤成一條線，然後在這條線上取某種均勻分配。

針對第二種方法，我們是「針對一個面取隨機數」，比起第一種方式，可以想成我們對一個 multivariate uniform distribution 取隨機數，所以我們是把圓面積透過某種映射方式映到二維的均勻分配。

針對第三種方式，我們是在一條半徑上取隨機數，與第一種一維的均勻分配比較相像，但結果卻大不相同。

這也是 3Blue1Brown 頻道作者另外強調的，我們看見的這個悖論其實是來自「我們不知道一條弦真正的定義是什麼」，所以不同的映射方式（上面每一種方法就對應到一種把圓映射到均勻分配的方式），就對應的不同的答案。

在影片中， numberphile 頻道的作者問道：「如果把所有可能的弦都蒐集起來，並一一編號，這樣是不是就可以解決 ambiguity 的問題？」然而，我們曉得所有的弦是不可數無窮多的，所以至少沒辦法用自然數幫他們標號。當我們使用實數（或者 [0,1] ）幫他們編號時，我們就會遇到「要怎麼編」的問題，而這個過程就會牽涉到我們具體而言把我們的「弦的分布」映到怎麼樣的分配裡去了。

R Code

# You have to be careful when the probabilities with an infinite space are involved
# What is the probability of the average length of the chord (say l) decided greater than an equilateral triangle?
# Prob( l > SOME_NUMBER) = ?

library(plotrix)
# Define an unite circle
plotInit <- function(){
  plot(x = c(1, 0, -1, 0), 
       y = c(0, 1, 0, -1),
       asp = 1, type = 'n', 
       xlim = c(-1, 1), ylim = c(-1, 1))
  draw.circle(0, 0, 1)
  lines(c(0, sqrt(3)/2, -sqrt(3)/2, 0), c(1, -1/2, -1/2, 1), col = "red")
}
# define transparency of lines
alpha <- rgb(red = 0, green = 0, blue = 1, alpha = 0.1)


# Perspective 1
# random sampling on the circle (在圓上隨機任取兩點，問通過此二點的弦長大於圓內接正三角形邊長的機率為何)
plotInit()
chord_len <- c()
for(i in 1:2000){
  theta <- runif(2, 0, 4*pi)
  X0 <- 1*cos(theta[1])
  Y0 <- 1*sin(theta[1])
  X1 <- 1*cos(theta[2])
  Y1 <- 1*sin(theta[2])
  chord_len <- c(chord_len, sqrt((X0-X1)^2 + (Y0-Y1)^2))
  # points(X0, Y0, col = alpha)
  # points(X1, Y1, col = alpha)
  lines(c(X0, X1), c(Y0, Y1), col = alpha)
}
sum(chord_len > sqrt(3))/length(chord_len)


# Perspective 2
# random sampling on the area of the circle
plotInit()
chord_len <- c()
for(i in 1:1000){
  X0 <- runif(1, -1, 1)
  Y0 <- runif(1, -1, 1)
  while( X0^2 + Y0^2 >= 1 ){
    X0 <- runif(1, -1, 1)
    Y0 <- runif(1, -1, 1)
  }
  points(X0, Y0, col = alpha)
  # lines(c(0, X0), c(0, Y0), col = "blue") # 圓心到圓內部一點
  k <- (1-X0^2-Y0^2)^(1/2)/sqrt(X0^2+Y0^2) # scale factor of normal vector
  # lines(c(X0, X0-k*Y0), c(Y0, Y0+k*X0), col = alpha) # 找出與此線段垂直的弦的一半
  # lines(c(X0, X0+k*Y0), c(Y0, Y0-k*X0), col = alpha) # 找出與此線段垂直的弦的另一半
  lines(c(X0-k*Y0, X0+k*Y0), c(Y0+k*X0, Y0-k*X0), col = alpha)
  # 弦之兩端點座標為 a(X0-k*Y0, Y0+k*X0); b(X0+k*Y0, Y0-k*X0)
  chord_len <- c(chord_len, sqrt(((X0-k*Y0) - (X0+k*Y0))^2 + ((Y0+k*X0) - (Y0-k*X0))^2))
}
sum(chord_len > sqrt(3))/length(chord_len)


# Perspective 3
# random sampling on the chord
plotInit()
chord_len <- c()
for(i in 1:1000){
  theta <- runif(1, 0, 4*pi)
  X0 <- 1*cos(theta)
  Y0 <- 1*sin(theta)
  # points(X0, Y0, col = "blue")
  # lines(c(0, X0), c(0, Y0), col = "blue") # 圓心到弦
  # 在弦上隨機選一點
  k <- runif(1, 0, sqrt(X0^2+Y0^2)) # 圓心到弦的線段長度為一，在此線段上抽
  
  X1 <- k*X0
  Y1 <- k*Y0
  # points(X1, Y1, col = "red") # 抽出一點，找過此點的弦，並垂直於此點至圓心的連線
  kprime <- (1-X1^2-Y1^2)^(1/2)/sqrt(X1^2+Y1^2) # scale factor of normal vector
  lines(c(X1, X1-kprime*Y1), c(Y1, Y1+kprime*X1), col = alpha) # 找出與此線段垂直的弦的一半
  lines(c(X1, X1+kprime*Y1), c(Y1, Y1-kprime*X1), col = alpha) # 找出與此線段垂直的弦的另一半
  # 弦之兩端點座標為 a(X1-k*Y1, Y1+k*X1); b(X1+k*Y1, Y1-k*X1)
  chord_len <- c(chord_len, sqrt(((X1-kprime*Y1) - (X1+kprime*Y1))^2 + ((Y1+kprime*X1) - (Y1-kprime*X1))^2))
}
sum(chord_len > sqrt(3))/length(chord_len)