伯特兰悖论模拟
Bertrand's Paradox Simulation
所以，我前几天在 Hacker News 上看到了这个：http://web.mit.edu/tee/www/bertrand/problem.html
它基本上说的是半径为 1 的圆上的随机弦的长度大于 3 的平方根的概率。
看起来答案是1/3似乎很明显，但是HN上的评论有比我聪明的人在争论这个问题。 https://news.ycombinator.com/item?id=10000926
我不想辩论，但我确实想确保我没有疯。所以我编码了我认为可以证明它是 P = 1/3 的代码，但我最终得到 P ~ .36。所以，我的代码一定有问题。
我可以进行完整性检查吗？
package com.jonas.betrand;
import java.awt.geom.Point2D;
import java.util.Random;
public class Paradox {
final static double ROOT_THREE = Math.sqrt(3);
public static void main(String[] args) {
int greater = 0;
int less = 0;
for (int i = 0; i < 1000000; i++) {
Point2D.Double a = getRandomPoint();
Point2D.Double b = getRandomPoint();
//pythagorean
if (Math.sqrt(Math.pow((a.x - b.x), 2) + Math.pow((a.y - b.y), 2)) > ROOT_THREE) {
greater++;
} else {
less++;
}
}
System.out.println("Probability Observerd: " + (double)greater/(greater+less));
}
public static Point2D.Double getRandomPoint() {
//get an x such that -1 < x < 1
double x = Math.random();
boolean xsign = new Random().nextBoolean();
if (!xsign) {
x *= -1;
}
//formula for a circle centered on origin with radius 1: x^2 + y^2 = 1
double y = Math.sqrt(1 - (Math.pow(x, 2)));
boolean ysign = new Random().nextBoolean();
if (!ysign) {
y *= -1;
}
Point2D.Double point = new Point2D.Double(x, y);
return point;
}
}
编辑：感谢一群人让我直截了当，我发现我找到随机点的方法确实不是那么随机。这是 returns 大约 1/3 的函数修复。
public static Point2D.Double getRandomPoint() {
//get an x such that -1 < x < 1
double x = Math.random();
Random r = new Random();
if (!r.nextBoolean()) {
x *= -1;
}
//circle centered on origin: x^2 + y^2 = r^2. r is 1.
double y = Math.sqrt(1 - (Math.pow(x, 2)));
if (!r.nextBoolean()) {
y *= -1;
}
if (r.nextBoolean()) {
return new Point2D.Double(x, y);
} else {
return new Point2D.Double(y, x);
}
}
伯特兰悖论就是这样：一个悖论。根据对问题的解释方式，答案可能是 1/3 或 1/2。看起来你采用了随机和弦方法，其中线的一侧是固定的，然后你将随机和弦绘制到圆的任何部分。使用这种方法，画出比sqrt(3)长的和弦的几率确实是1/3.
但是如果你使用不同的方法，我会称之为随机半径方法，你会发现它可以是 1/2！随机半径是这样的，你在圆圈里画一个半径，然后你取这个半径平分的随机弦。此时随机一个和弦会比sqrt(3)长1/2的时候。
最后，随机中点法。在圆圈中随机选择一个点，然后以这个随机点作为弦的中点绘制弦。如果该点落在半径为 1/2 的同心圆内，则弦比 sqrt(3) 短。如果落在同心圆外，则比sqrt(3)长。半径为 1/2 的圆的面积是半径为 1 的圆的 1/4，所以弦小于 sqrt(3) 的概率是 1/4.
至于你的代码，我还没来得及看，但希望这能澄清悖论（这只是一个不完整的问题，实际上不是悖论）:D
我相信你需要假设一个固定点在 (0, 1) .
真见鬼，我在 Swift 中写了错误的版本（学习 Swift！）：
struct P {
let x, y: Double
init() {
x = (Double(arc4random()) / 0xFFFFFFFF) * 2 - 1
y = sqrt(1 - x * x) * (arc4random() % 2 == 0 ? 1 : -1)
}
func dist(other: P) -> Double {
return sqrt((x - other.x) * (x - other.x) + (y - other.y) * (y - other.y))
}
}
let root3 = sqrt(3.0)
let total = 100_000_000
var samples = 0
for var i = 0; i < total; i++ {
if P().dist(P()) > root3 {
samples++
}
}
println(Double(samples) / Double(total))
答案确实是0.36。正如评论所解释的那样，随机 X 值更有可能选择 pi/2 附近的 "flattened area" 并且极不可能选择 0 和 pi 附近的 "vertically squeezed" 区域。
但是在 P 的构造函数中很容易修复：
(Double(arc4random()) / 0xFFFFFFFF 是 [0, 1) 中的随机浮点数的花哨说法)
let angle = Double(arc4random()) / 0xFFFFFFFF * M_PI * 2
x = cos(angle)
y = sin(angle)
// outputs 0.33334509
我认为伯特兰悖论与其说是悖论，不如说是概率论中的警示。真的是在问这个问题：random是什么意思？
Bertrand 认为随机选择和弦有三种自然但不同的方法，给出了三个不同的答案。但当然，还有其他随机方法，但这些方法可以说不是最自然的方法（也就是说，不是第一个想到的）。例如，我们可以以 non-uniform 的方式随机放置两个弦端点。或者我们根据一些 non-uniform 密度定位和弦中点，例如截断的 bi-variate 法线。
要用编程语言模拟这三种方法，需要能够在单位区间内生成统一的随机变量，这是所有标准（伪）随机数生成器都应该做的。对于 methods/solutions 之一（随机中点之一），您必须取其中一个均匀随机变量的平方根。然后，您将随机变量乘以一个合适的因子（或重新缩放）。然后对于每种模拟方法（或解决方案），一些几何给出了两个端点的表达式。
有关更多详细信息，我已经编写了一套 post about this problem. I recommend the links and books I have cited at the end of that post, under the section Further reading. For example, see Section 1.3 in this 新的已发布讲义。伯特兰悖论也出现在 Isaac 的 The Pleasures of Probability 中。 Clark 在 Paradoxes from A to Z 一书中以 non-mathematical 的方式进行了介绍。
我也上传了一些MATLAB、R和Python中的仿真代码，可以找到here。
例如，在 Python 中（使用 NumPy）：
import numpy as np; #NumPy package for arrays, random number generation, etc
import matplotlib.pyplot as plt #for plotting
from matplotlib import collections as mc #for plotting line chords
###START Parameters START###
#Simulation disk dimensions
xx0=0; yy0=0; #center of disk
r=1; #disk radius
numbLines=10**2;#number of lines
###END Parameters END###
###START Simulate three solutions on a disk START###
#Solution A
thetaA1=2*np.pi*np.random.uniform(0,1,numbLines); #choose angular component uniformly
thetaA2=2*np.pi*np.random.uniform(0,1,numbLines); #choose angular component uniformly
#calculate chord endpoints
xxA1=xx0+r*np.cos(thetaA1);
yyA1=yy0+r*np.sin(thetaA1);
xxA2=xx0+r*np.cos(thetaA2);
yyA2=yy0+r*np.sin(thetaA2);
#calculate midpoints of chords
xxA0=(xxA1+xxA2)/2; yyA0=(yyA1+yyA2)/2;
#Solution B
thetaB=2*np.pi*np.random.uniform(0,1,numbLines); #choose angular component uniformly
pB=r*np.random.uniform(0,1,numbLines); #choose radial component uniformly
qB=np.sqrt(r**2-pB**2); #distance to circle edge (alonge line)
#calculate trig values
sin_thetaB=np.sin(thetaB);
cos_thetaB=np.cos(thetaB);
#calculate chord endpoints
xxB1=xx0+pB*cos_thetaB+qB*sin_thetaB;
yyB1=yy0+pB*sin_thetaB-qB*cos_thetaB;
xxB2=xx0+pB*cos_thetaB-qB*sin_thetaB;
yyB2=yy0+pB*sin_thetaB+qB*cos_thetaB;
#calculate midpoints of chords
xxB0=(xxB1+xxB2)/2; yyB0=(yyB1+yyB2)/2;
#Solution C
#choose a point uniformly in the disk
thetaC=2*np.pi*np.random.uniform(0,1,numbLines); #choose angular component uniformly
pC=r*np.sqrt(np.random.uniform(0,1,numbLines)); #choose radial component
qC=np.sqrt(r**2-pC**2); #distance to circle edge (alonge line)
#calculate trig values
sin_thetaC=np.sin(thetaC);
cos_thetaC=np.cos(thetaC);
#calculate chord endpoints
xxC1=xx0+pC*cos_thetaC+qC*sin_thetaC;
yyC1=yy0+pC*sin_thetaC-qC*cos_thetaC;
xxC2=xx0+pC*cos_thetaC-qC*sin_thetaC;
yyC2=yy0+pC*sin_thetaC+qC*cos_thetaC;
#calculate midpoints of chords
xxC0=(xxC1+xxC2)/2; yyC0=(yyC1+yyC2)/2;
###END Simulate three solutions on a disk END###
所以，我前几天在 Hacker News 上看到了这个：http://web.mit.edu/tee/www/bertrand/problem.html
它基本上说的是半径为 1 的圆上的随机弦的长度大于 3 的平方根的概率。
看起来答案是1/3似乎很明显，但是HN上的评论有比我聪明的人在争论这个问题。 https://news.ycombinator.com/item?id=10000926
我不想辩论，但我确实想确保我没有疯。所以我编码了我认为可以证明它是 P = 1/3 的代码，但我最终得到 P ~ .36。所以，我的代码一定有问题。
我可以进行完整性检查吗？
package com.jonas.betrand;
import java.awt.geom.Point2D;
import java.util.Random;
public class Paradox {
final static double ROOT_THREE = Math.sqrt(3);
public static void main(String[] args) {
int greater = 0;
int less = 0;
for (int i = 0; i < 1000000; i++) {
Point2D.Double a = getRandomPoint();
Point2D.Double b = getRandomPoint();
//pythagorean
if (Math.sqrt(Math.pow((a.x - b.x), 2) + Math.pow((a.y - b.y), 2)) > ROOT_THREE) {
greater++;
} else {
less++;
}
}
System.out.println("Probability Observerd: " + (double)greater/(greater+less));
}
public static Point2D.Double getRandomPoint() {
//get an x such that -1 < x < 1
double x = Math.random();
boolean xsign = new Random().nextBoolean();
if (!xsign) {
x *= -1;
}
//formula for a circle centered on origin with radius 1: x^2 + y^2 = 1
double y = Math.sqrt(1 - (Math.pow(x, 2)));
boolean ysign = new Random().nextBoolean();
if (!ysign) {
y *= -1;
}
Point2D.Double point = new Point2D.Double(x, y);
return point;
}
}
编辑：感谢一群人让我直截了当，我发现我找到随机点的方法确实不是那么随机。这是 returns 大约 1/3 的函数修复。
public static Point2D.Double getRandomPoint() {
//get an x such that -1 < x < 1
double x = Math.random();
Random r = new Random();
if (!r.nextBoolean()) {
x *= -1;
}
//circle centered on origin: x^2 + y^2 = r^2. r is 1.
double y = Math.sqrt(1 - (Math.pow(x, 2)));
if (!r.nextBoolean()) {
y *= -1;
}
if (r.nextBoolean()) {
return new Point2D.Double(x, y);
} else {
return new Point2D.Double(y, x);
}
}
伯特兰悖论就是这样：一个悖论。根据对问题的解释方式，答案可能是 1/3 或 1/2。看起来你采用了随机和弦方法，其中线的一侧是固定的，然后你将随机和弦绘制到圆的任何部分。使用这种方法，画出比sqrt(3)长的和弦的几率确实是1/3.
但是如果你使用不同的方法，我会称之为随机半径方法，你会发现它可以是 1/2！随机半径是这样的，你在圆圈里画一个半径，然后你取这个半径平分的随机弦。此时随机一个和弦会比sqrt(3)长1/2的时候。
最后，随机中点法。在圆圈中随机选择一个点，然后以这个随机点作为弦的中点绘制弦。如果该点落在半径为 1/2 的同心圆内，则弦比 sqrt(3) 短。如果落在同心圆外，则比sqrt(3)长。半径为 1/2 的圆的面积是半径为 1 的圆的 1/4，所以弦小于 sqrt(3) 的概率是 1/4.
至于你的代码，我还没来得及看，但希望这能澄清悖论（这只是一个不完整的问题，实际上不是悖论）:D
我相信你需要假设一个固定点在 (0, 1) .
真见鬼，我在 Swift 中写了错误的版本（学习 Swift！）：
struct P {
let x, y: Double
init() {
x = (Double(arc4random()) / 0xFFFFFFFF) * 2 - 1
y = sqrt(1 - x * x) * (arc4random() % 2 == 0 ? 1 : -1)
}
func dist(other: P) -> Double {
return sqrt((x - other.x) * (x - other.x) + (y - other.y) * (y - other.y))
}
}
let root3 = sqrt(3.0)
let total = 100_000_000
var samples = 0
for var i = 0; i < total; i++ {
if P().dist(P()) > root3 {
samples++
}
}
println(Double(samples) / Double(total))
答案确实是0.36。正如评论所解释的那样，随机 X 值更有可能选择 pi/2 附近的 "flattened area" 并且极不可能选择 0 和 pi 附近的 "vertically squeezed" 区域。
但是在 P 的构造函数中很容易修复：
(Double(arc4random()) / 0xFFFFFFFF 是 [0, 1) 中的随机浮点数的花哨说法)
let angle = Double(arc4random()) / 0xFFFFFFFF * M_PI * 2
x = cos(angle)
y = sin(angle)
// outputs 0.33334509
我认为伯特兰悖论与其说是悖论，不如说是概率论中的警示。真的是在问这个问题：random是什么意思？
Bertrand 认为随机选择和弦有三种自然但不同的方法，给出了三个不同的答案。但当然，还有其他随机方法，但这些方法可以说不是最自然的方法（也就是说，不是第一个想到的）。例如，我们可以以 non-uniform 的方式随机放置两个弦端点。或者我们根据一些 non-uniform 密度定位和弦中点，例如截断的 bi-variate 法线。
要用编程语言模拟这三种方法，需要能够在单位区间内生成统一的随机变量，这是所有标准（伪）随机数生成器都应该做的。对于 methods/solutions 之一（随机中点之一），您必须取其中一个均匀随机变量的平方根。然后，您将随机变量乘以一个合适的因子（或重新缩放）。然后对于每种模拟方法（或解决方案），一些几何给出了两个端点的表达式。
有关更多详细信息，我已经编写了一套 post about this problem. I recommend the links and books I have cited at the end of that post, under the section Further reading. For example, see Section 1.3 in this 新的已发布讲义。伯特兰悖论也出现在 Isaac 的 The Pleasures of Probability 中。 Clark 在 Paradoxes from A to Z 一书中以 non-mathematical 的方式进行了介绍。
我也上传了一些MATLAB、R和Python中的仿真代码，可以找到here。
例如，在 Python 中（使用 NumPy）：
import numpy as np; #NumPy package for arrays, random number generation, etc
import matplotlib.pyplot as plt #for plotting
from matplotlib import collections as mc #for plotting line chords
###START Parameters START###
#Simulation disk dimensions
xx0=0; yy0=0; #center of disk
r=1; #disk radius
numbLines=10**2;#number of lines
###END Parameters END###
###START Simulate three solutions on a disk START###
#Solution A
thetaA1=2*np.pi*np.random.uniform(0,1,numbLines); #choose angular component uniformly
thetaA2=2*np.pi*np.random.uniform(0,1,numbLines); #choose angular component uniformly
#calculate chord endpoints
xxA1=xx0+r*np.cos(thetaA1);
yyA1=yy0+r*np.sin(thetaA1);
xxA2=xx0+r*np.cos(thetaA2);
yyA2=yy0+r*np.sin(thetaA2);
#calculate midpoints of chords
xxA0=(xxA1+xxA2)/2; yyA0=(yyA1+yyA2)/2;
#Solution B
thetaB=2*np.pi*np.random.uniform(0,1,numbLines); #choose angular component uniformly
pB=r*np.random.uniform(0,1,numbLines); #choose radial component uniformly
qB=np.sqrt(r**2-pB**2); #distance to circle edge (alonge line)
#calculate trig values
sin_thetaB=np.sin(thetaB);
cos_thetaB=np.cos(thetaB);
#calculate chord endpoints
xxB1=xx0+pB*cos_thetaB+qB*sin_thetaB;
yyB1=yy0+pB*sin_thetaB-qB*cos_thetaB;
xxB2=xx0+pB*cos_thetaB-qB*sin_thetaB;
yyB2=yy0+pB*sin_thetaB+qB*cos_thetaB;
#calculate midpoints of chords
xxB0=(xxB1+xxB2)/2; yyB0=(yyB1+yyB2)/2;
#Solution C
#choose a point uniformly in the disk
thetaC=2*np.pi*np.random.uniform(0,1,numbLines); #choose angular component uniformly
pC=r*np.sqrt(np.random.uniform(0,1,numbLines)); #choose radial component
qC=np.sqrt(r**2-pC**2); #distance to circle edge (alonge line)
#calculate trig values
sin_thetaC=np.sin(thetaC);
cos_thetaC=np.cos(thetaC);
#calculate chord endpoints
xxC1=xx0+pC*cos_thetaC+qC*sin_thetaC;
yyC1=yy0+pC*sin_thetaC-qC*cos_thetaC;
xxC2=xx0+pC*cos_thetaC-qC*sin_thetaC;
yyC2=yy0+pC*sin_thetaC+qC*cos_thetaC;
#calculate midpoints of chords
xxC0=(xxC1+xxC2)/2; yyC0=(yyC1+yyC2)/2;
###END Simulate three solutions on a disk END###