4.4 Bertrand’s Paradox in Probability

Joseph Bertrand라는 프랑스 수학자가 1889년에 제시한 문제이다. 자세한 내력은 https://en.wikipedia.org/wiki/Bertrand_paradox_(probability)를 참고한다.

평면상에 원이 있고, 이에 내접하는 정삼각형이 있다고 하자. 이때 직선을 임의로 그었을 때, 할선(chord, 직선과 원의 교점을 잇는 선분)의 길이가 정삼각형 한 변의 길이보다 클 확률을 구하라는 문제이다.

얼핏 보면 문제가 잘 정의되어 있어서, 답이 명확한 것 같지만, 여기에는 3가지의 똑같이 믿을 만한 답(1/2, 1/3, 1/4)이 있다.

방법 1: 원의 둘레 상의 임의의 두 점을 잡아 선분을 만든다. 이 때 두 점과 원의 중심을 연결한 선분들이 원의 중심에서 이루는 각도가 (반시계방향으로 회전하는 각도로 생각하면 0에서 360도까지의 범위이고 이 때) 120에서 240도일 때이다. (작은 각만을 기준으로 할 때는 0에서 180도까지의 범위이고 이 때) 120에서 180도까지 일 때이다. 따라서, 확률은 (240 - 120)/360 또는 (180 - 120)/180, 즉, 1/3이다.

방법 2: 원의 중심과 선분의 중심의 거리가 원의 반지름의 반보다 가까울 때이므로 1/2이다.

방법 3: 선분의 중심이 원래의 원보다 반지름이 반이고, 중심이 같은 원 안에 놓일 때이므로 면적비로 하여 1/4이다.

각각에 대하여 simulation을 해보면 다음과 같다.

4.4.1 방법 1: 원 둘레 상의 임의의 두 점을 잡을 때

N = 2000
Radius = 1
th = seq(0, 2*pi, length.out=(N + 1))

DistB = matrix(NA, nrow=N, ncol=3)
for (i in 1:N) {
  DistB[i, 1:2] = sample(th, 2, replace=T)
}

DistB[, 3] = abs(DistB[,1] - DistB[,2])
mean(DistB[, 3] >= 2/3*pi & DistB[, 3] <= 4/3*pi) # 1/3

[1] 0.3445

### Chords plot
#dev.new()
oPar = par(list(mfrow=c(1, 2)))
plot(0, 0, type="n", xlim=c(-1, 1), ylim=c(-1, 1), asp=1, bty="n", xlab="Chords",
     ylab="", xaxp=c(-1, 1, 2), yaxp=c(-1, 1, 2))
for (i in 1:200) {
  if (DistB[i, 3] >= 2/3*pi & DistB[i, 3] <= 4/3*pi) {
    lines(cos(DistB[i, 1:2]), sin(DistB[i, 1:2]), col="red")
  } else {
    lines(cos(DistB[i, 1:2]), sin(DistB[i, 1:2]), col="blue")
  }
}
lines(cos(th), sin(th), lwd=2)

### Midpoints plot
#dev.new()
midXs = Radius*(cos(DistB[, 1]) + cos(DistB[, 2]))/2
midYs = Radius*(sin(DistB[, 1]) + sin(DistB[, 2]))/2
plot(midXs, midYs, type="p", xlim=c(-1, 1), ylim=c(-1, 1), asp=1, pch=20, cex=0.1,
     bty="n", yaxt="n", xlab="Midpoints", ylab="", xaxp=c(-1, 1, 2), yaxp=c(-1, 1, 2))
lines(cos(th), sin(th), lwd=2)

Figure 4.3: Method 1 - Two random endpoints

par(oPar)

4.4.2 방법 2: 원의 중심과 선분의 중점 사이 거리

아래에서 Min과 Max의 절대값을 크게 할 수록 더 정확하지만 시간이 훨씬 오래 걸리게 된다.

Max = 100
Min = -100

Dist = matrix(NA, nrow=N, ncol=4)
n = 0
while (n < N) {
  x1 = runif(1, Min, Max)
  y1 = runif(1, Min, Max)
  m = tan(runif(1, 0, pi))
  dist1 = abs(m*x1 - y1)/sqrt(m*m + 1)
  if (dist1 > Radius) next
  n = n + 1
  Dist[n, ] = c(x1, y1, m, dist1)
}

mean(Dist[, 4] < Radius/2) # 0.5

[1] 0.4955

### Chords plot
#dev.new()
oPar = par(list(mfrow=c(1, 2)))
plot(0, 0, type="n", xlim=c(-1, 1), ylim=c(-1, 1), asp=1, bty="n", xlab="Chords",
     ylab="", xaxp=c(-1, 1, 2), yaxp=c(-1, 1, 2))
for (i in 1:200) {
  if (Dist[i, 4] <= Radius/2) {
    abline(a = -Dist[i, 1]*Dist[i, 3] + Dist[i, 2], b = Dist[i, 3], col="red")
  } else {
    abline(a = -Dist[i, 1]*Dist[i, 3] + Dist[i, 2], b = Dist[i, 3], col="blue")
  }
}
lines(cos(th), sin(th), lwd=2)

### Midpoints plot
ms = Dist[, 3]
ns = -Dist[, 3]*Dist[, 1] + Dist[, 2]
msqp1 = ms*ms + 1 # m square + 1
Xs = -ms*ns/msqp1
Ys = ns/msqp1
#dev.new()
plot(Xs, Ys, type="p", xlim=c(-1, 1), ylim=c(-1, 1), asp=1, pch=20, cex=0.1, bty="n",
     yaxt="n", xlab="Midpoints", ylab="", xaxp=c(-1, 1, 2), yaxp=c(-1, 1, 2))
lines(cos(th), sin(th), lwd=2)

Figure 4.4: Method 2 - Random x1, y1, and theta

par(oPar)

4.4.3 방법 3: 선분의 중점이 (원래 원의 중심과 같고) 반지름이 반인 원 안에 있을 때

DistC = matrix(NA, nrow=N, ncol=3)
for (i in 1:N) {
  DistC[i, 1:2] = runif(2, -1, 1)
  DistC[i, 3] = sqrt(DistC[i, 1]^2 + DistC[i, 2]^2)
}

Dist1 = DistC[DistC[,3] < 1,]
mean(Dist1[,3] < Radius/2) # 1/4

[1] 0.2438252

### Chords plot
ms = tan(runif(N, 0, pi))
ns = -ms*DistC[, 1] + DistC[, 2]
#dev.new()
oPar = par(list(mfrow=c(1, 2)))
plot(0, 0, type="n", xlim=c(-1, 1), ylim=c(-1, 1), asp=1, bty="n", xlab="Chords",
     ylab="", xaxp=c(-1, 1, 2), yaxp=c(-1, 1, 2))
for (i in 1:200) {
  if (DistC[i, 3] <= Radius/2) {
    abline(a = ns[i], b = ms[i], col="red") # Note this is Dist not DistC
  } else {
    abline(a = ns[i], b = ms[i], col="blue") # Note this is Dist not DistC
  }
}
lines(cos(th), sin(th), lwd=2)

### Midpoints plot
#dev.new()
plot(DistC[, 1], DistC[, 2], type="p", xlim=c(-1, 1), ylim=c(-1, 1), asp=1, pch=20, bty="n",
     cex=0.1, yaxt="n", xlab="Midpoints", ylab="", xaxp=c(-1, 1, 2), yaxp=c(-1, 1, 2))
lines(cos(th), sin(th), lwd=2)
lines(cos(th)/2, sin(th)/2, lwd=2)

Figure 4.5: Method 3 - Random midpoint

par(oPar)

이 문제를 다시 정리하자면 방법 1은 극좌표계(polar coordinate)를 쓴 것이며, 방법 2는 1차원 좌표계, 방법 3은 2차원 직교좌표계(Cartesian coordinate)를 쓴 것이다. 일반적으로 1차원 좌표계를 쓴 경우의 제곱(1/2의 제곱)이 2차원 직교좌표계의 답이 된다. 여기에서의 교훈은 어떤 좌표계를 쓰냐에 따라 답이 다를 수 있다는 점에 유의하라는 것이다.

Slope = 1
mean(tan(runif(10000, 0, pi/2)) < Slope) # Theoretically 1 - atan(0.5)/(pi/4) = 0.4096655

[1] 0.4896

atan(Slope) # 0.7853982 = pi/4

[1] 0.7853982

atan(Slope)/(pi/2) # 0.5

[1] 0.5

mp = tan(pi/4) # 0.4142136
mean(tan(runif(10000, 0, pi/2)) > mp) # 0.5

[1] 0.4975

plot(0, 0, type="n", xlim=c(-1, 1), ylim=c(-1, 1), asp=1, xlab="", ylab="")
for (i in 1:9) {
  abline(a = 0, b = i*1e6, col="red")
  abline(a = 0, b = tan(i/10*pi/2), col="blue")
}
abline(h = 0)
abline(a = 0, b = 1)
abline(v = c(-1, 1))

y0 = seq(0, 1, length.out=6)
points(x = c(rep(1, 6), rep(-1, 6)), y = c(y0, -y0), pch=20, col="red")
points(x = cos((0:10)/10*pi/2), y = sin((0:10)/10*pi/2), pch=20, col="blue")
th = seq(0, 2*pi, length.out=101)
lines(cos(th), sin(th), lwd=2)

Figure 4.6: Even thetas vs. even y coordinates

직교좌표계이냐, 극좌표계이냐에 따라 점들의 밀도(특정 점이 골라질 확률)가 다르게 분포됨을 알 수 있다.

plot(0, 0, type="n", xlim=c(-1, 1), ylim=c(-1, 1), asp=1, xlab="", ylab="")

rs = seq(0, 1, length.out=11)
th = seq(0, 2*pi, length.out=41)
xs = vector(length=11*41)
ys = xs
for (i in 1:(11*41)) {
  j = (i %/% 41) + 1
  k = (i %% 41) + 1
  xs[i] = rs[j]*cos(th[k])
  ys[i] = rs[j]*sin(th[k])
  abline(a = 0, b = tan(th[k]), col="red")
}
points(xs, ys, pch=20, col="red")

xs = rep(seq(-1.2, 1.2, length.out=25), 25)
ys = sort(rep(seq(-1.2, 1.2, length.out=25), 25))
points(xs, ys, pch=20)

th = seq(0, 2*pi, length.out=101)
for (r in (1:15)/10) {
  lines(r*cos(th), r*sin(th), lwd=2, col="red")
}

Figure 4.7: Polar coordinates vs. Cartesian coordinates