counterintuitive359jaynes.py

# Jaynes' counter-intuitive equation
# 3.59 in PT:TLoS

import random

random.seed(1) # For experiment reproduction

# P(R1|R2)
count_reds1 = 0.0
count_total1 = 0.0
p1 = []
p1exp = []
# P(R1|RLater)
count_reds2 = 0.0
count_total2 = 0.0
p2 = []
p2exp = []

# Repeat the experiment a lot of times
experiments = []
exp = []
for i in range(1000000): # Experiment repetition
    experiments.append([])
    # Four balls, two of them red
    n_red = 2.0
    n_white = 2.0
    for j in range(3): # Three draws
        if random.random() < n_red/(n_red + n_white):
            n_red -= 1.0
            experiments[i].append(True)
        else:
            n_white -= 1.0
            experiments[i].append(False)
    urn = experiments[i]
    if urn[1]: # Second draw is red
        count_reds1 += urn[0] # Count times first draw is red
        count_total1 += 1.0
        p1.append(count_reds1/count_total1)
        p1exp.append(i)
    if urn[1] or urn[2]: # Second or third draw is red
        count_reds2 += urn[0] # Count times first draw is red
        count_total2 += 1.0
        p2.append(count_reds2/count_total2)
        p2exp.append(i)

# Knowing that the second draw is red, what is the probability of the first being red? Jaynes' answer: 1/3
print("P(R1|R2) = " + str(count_reds1/count_total1) + " with size " + str(len(p1)))
# Knowing that there will be at least one red in the next two draws, what is the probability the first one is a red? Jaynes' answer: 2/5
print("P(R1|RLater) = " + str(count_reds2/count_total2) + " with size " + str(len(p2)))

import pylab
pylab.plot(p1, p1exp)
pylab.plot(p2, p2exp)
pylab.show()