javasimulationparadox

Bertrand's Paradox Simulation


So, I saw this on Hacker News the other day: http://web.mit.edu/tee/www/bertrand/problem.html

It basically says what's the probability that a random chord on a circle with radius of 1 has a length greater than the square root of 3.

Looking at it, it seems obvious that the answer is 1/3, but comments on HN have people who are smarter than me debating this. https://news.ycombinator.com/item?id=10000926

I didn't want to debate, but I did want to make sure I wasn't crazy. So I coded what I thought would prove it to be P = 1/3, but I end up getting P ~ .36. So, something's got to be 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 bunch of people setting me straight, I found that my method of finding a random point wasn't indeed so random. Here is a fix for that 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);
            }
        }

Solution

  • I believe you need to assume one fixed point say at (0, 1) and then choose a random amount of rotation in [0, 2*pi] around the circle for the location of the second point of the chord.

    Just for the hell of it I wrote your incorrect 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 the comments have been explaining, a random X value is more likely to choose the "flattened area" around pi/2 and highly unlikely to choose the "vertically squeezed" area around 0 and pi.

    It is easily fixed however in the constructor for P: (Double(arc4random()) / 0xFFFFFFFF is fancy-speak for random floating point number in [0, 1))

            let angle = Double(arc4random()) / 0xFFFFFFFF * M_PI * 2
            x = cos(angle)
            y = sin(angle)
            // outputs 0.33334509