birthday problem for 4 people

bd.py

import numpy
import matplotlib.pyplot as plt

def prob(M):
	# for M people, compute the probability of having more than 4 with same birthday
	hits = 0
	# number of simulation instances
	N = 1000
	I = numpy.arange(365).reshape((1,-1))
	for j in range(N):
		s = numpy.random.randint(0, 365, size=M)
		hit2 = ((s.reshape((-1,1)) == I).sum(axis=0) >= 4).any()
		hits += hit2
	return hits * 1. / N

A = 1
O = 300

# bisection
while True:
	M = (A + O) / 2
	if M == A or M == O: 
		break
	p = prob(M)
	print '%4d %.2f%%' % (M, p*100)
	plt.plot(M, p, 'x')
	plt.savefig('bd.pdf', bbox_inches='tight')
	
	if p > 0.5:
		O = M
	elif p < 0.5:
		A = M

print 'solution is between %d and %d' % (A, O)

x = numpy.arange(100, 300, 10)
y = [prob(M) for M in x]
plt.plot(x, y)
plt.savefig('bd.pdf', bbox_inches='tight')
plt.savefig('bd.png', bbox_inches='tight')

# Output: 
"""
 $ python bd.py 
 150 25.90%
 225 76.10%
 187 50.00%
 187 49.60%
 206 64.50%
 196 58.30%
 191 52.30%
 189 51.30%
 188 49.60%
solution is between 188 and 189
 $ python bd.py 
 150 25.80%
 225 76.70%
 187 50.00%
 187 50.70%
 168 37.60%
 177 45.00%
 182 47.30%
 184 46.20%
 185 46.70%
 186 50.40%
solution is between 185 and 186
"""