Modeling Bertrand Paradox
So, I saw this on Hacker News the other day: http://web.mit.edu/tee/www/bertrand/problem.html
Basically says that the probability that a random chord on a circle with radius 1 is longer than the square root of 3.
Looking at this, it seems obvious that the answer is 1/3, but there are people in the comments on HN who are smarter than I am discussing this. https://news.ycombinator.com/item?id=10000926
I didn't want to discuss, but I really wanted to make sure that I wasn't crazy. So I coded what I thought would prove it to be P = 1/3, but I get P ~ .36. So, something is wrong with my code.
Can I get a sanity check?
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;
}
}
EDIT: Thanks to a lot of people hiding me, I found that my method of finding a random point was not really that random. Here's a fix for this function, which returns about 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);
}
}
source to share
I suppose you need to take one fixed point, say at (0, 1), and then choose an arbitrary amount of rotation at [0, 2 * pi] around the circle for the location of the second chord point.
Just heck I wrote your wrong version in Swift (learn 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))
And the answer is indeed 0.36. As explained in the comments, a random value of X is more likely to choose a "flattened region" around pi / 2 and less likely to choose a "vertically collapsed" region around 0 and pi.
Easily fixed, however, in the constructor for P: ( Double(arc4random()) / 0xFFFFFFFF
is weird for a random floating point number in [0, 1))
let angle = Double(arc4random()) / 0xFFFFFFFF * M_PI * 2
x = cos(angle)
y = sin(angle)
// outputs 0.33334509
source to share
Bertrand's paradox is exactly this: a paradox. The answer can be asserted as 1/3 or 1/2 depending on how the problem is interpreted. You seem to pick a random chord where one side of the line is fixed, and then you produce a random chord anywhere in the circle. Using this method, the chances of drawing a chord that is longer than sqrt (3) are indeed 1/3.
But if you use a different approach, I will call it random radius, you will see that it could be 1/2! A random radius is this, you draw a radius in a circle and then take a random chord that divides that radius in half. At this point, the random chord will be longer than sqrt (3) 1/2 time.
Finally, the random middle method. Pick a random point in the circle and then draw a chord with that random point in the middle of the chord. If this point falls within a concentric circle of radius 1/2, then the chord is shorter than sqrt (3). If it falls outside the concentric circle, it is longer than sqrt (3). A circle of radius 1/2 has 1/4 the area of a circle with a radius of 1, so the probability that the chord is less than sqrt (3) is 1/4.
Regarding your code, I haven't had time to look at it yet, but hope this clears up the paradox (which is just an incomplete question, not actually a paradox): D
source to share
I would say that Bertrand's paradox is not a paradox, but rather a warning lesson in probability. The question is really asked here: what do you mean by random?
Bertrand argued that there are three natural but different methods for randomly choosing a chord, giving three different answers. But, of course, there are other random methods, but these methods are perhaps not the most natural (that is, not the first that come to mind). For example, we could arbitrarily arrange the two end points of the chords in an uneven way. Or we place the midpoint of the chord according to some non-uniform density, such as a truncated bidirectional normal.
To simulate the three methods in a programming language, you must be able to generate homogeneous random values on a unit interval, which is what all standard (pseudo) -random number generators should do. For one of the methods / solutions (random midpoint) you need to take the square root of one of the uniform random variables. Then you multiply the random values by a suitable factor (or change the scale). Then, for each modeling method (or solution), some geometry gives expressions for the two endpoints.
For more details, I wrote a post about this issue. I recommend the links and books I quoted at the end of this post in the Further Reading section. For example, see Section 1.3 in this new set of published lecture notes. Bertrand's paradox is also found in Isaac's Pleasures of Probability. This is non-mathematically highlighted in Clark's book "Paradoxes from A to Z".
I also downloaded the simulation code in MATLAB, R and Python which can be found here .
For example, in Python (with 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###
source to share