Python code to illustrate the Tuesday Boy paradox ("I have two children. One is a boy born on a Tuesday. What is the probability I have two boys") for a blog post.

paradox.py

# paradox.py
#
# Python code to illustrate the Tuesday Boy paradox ("I have two children. 
# One is a boy born on a Tuesday. What is the probability I have two boys") 
# for a blog post.
# 
# We're simulating two sampling modes to illustrate how the approach to 
# sampling, and posing the initial question, is important. This has general
# implications for what we can reasonably infer from experiments where we
# did not design the experiment to answer a specific question.
#
# (c) L.Pritchard 2013

import random                      # We need a uniform distribution
from collections import Counter    # Useful shortcut

# Define sexes and days
sexes = ['B', 'G']
days = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 
        'Saturday', 'Sunday']

# Fill a 'room' with 10000 parents, each having two children that are either
# boys (B) or girls (G), born on a specific day, indicated by a tuple
# (sex, day). This is a single realisation of uniform distributions of boys,
# girls, and days of birth
parents = []
for i in range(10000):
    children = ((random.choice(sexes), random.choice(days)),
                (random.choice(sexes), random.choice(days)))
    parents.append(children)

# Report the overall probabilities of having each sex
# This should give us P~0.25 for each category
categories = [(parent[0][0], parent[1][0]) for parent in parents]
counter = Counter(categories)
print '\n' + '\t'.join(["Category", "Count", "P"])
prob_dict = {}
for category in sorted(counter):
    prob_dict[category] = 1.*counter[category]/len(parents)
    print '\t'.join([str(e) for e in \
                         [category, counter[category], 
                          prob_dict[category]]])

print '\n' + '\t'.join(["Category", "Count", "P"])
for name, cats in [("Two girls:", [('G', 'G')]),
                   ("Two boys:", [('B', 'B')]),
                   ("One girl, one boy:", [('G', 'B'), ('B', 'G')])]:
    print '\t'.join([str(e) for e in [name, 
                                      sum([counter[c] for c in cats]),
                                      sum([prob_dict[c] for c in cats])]])

# Report the overall probabilities of having each combination of sex and 
# birthday
counter = Counter(parents)
print '\n' + '\t'.join(["Category", "Count", "P"])
for category in sorted(counter):
    prob_dict[category] = 1.*counter[category]/len(parents)
    print '\t'.join([str(e) for e in \
                         [category, counter[category], 
                          prob_dict[category]]])

# 1) Identify the proportion/probability of parents that satisfy the puzzle;
#    this involves stipulating that at least one boy in the final solution
#    set must have been born on a Tuesday
print "\nWe start with %d parents" % len(parents)
# Remove all parents that don't have at least one boy born on a Tuesday
soln1 = [p for p in parents if ('B', 'Tuesday') in p]
print "%d parents have at least one boy born on a Tuesday" % len(soln1)
soln2 = [p for p in soln1 if (p[0][0], p[1][0]) == ('B', 'B')]
print "Of these, %d have two boys" % len(soln2)
print "P(two boys|at least one boy born on a Tuesday) = %.3f" % \
    (1.*len(soln2)/len(soln1))

# 2) Identify the proportion/probability of parents that have two boys,
#    only after identifying the day on which a single boy in the sample
#    happens to have been born.
print "\nWe start with %d parents" % len(parents)
# Remove all parents that don't have at least one boy
soln1 = [p for p in parents if (p[0][0], p[1][0]) != ('G', 'G')]
print "%d parents have at least one boy" % len(soln1)
# Select someone to talk to:
parent = random.choice(soln1)
print "- And what day was (one of) your sons born on?"
for child in parent:
    if child[0] == 'B':
        day = child[1]
        print '- %s' % day
        break               # We only want the first one
soln2 = [p for p in soln1 if (p[0][0], p[1][0]) == ('B', 'B')]
print "Of the parents that are left, %d have two boys" % len(soln2)
print "P(two boys|at least one boy born on a %s) = %.3f" % \
    (day, 1.*len(soln2)/len(soln1))