#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

set.seed(0)

# Exercise 30:
#
birthday_problem = function(nSims=10,nPeopleInRoom=25,nPeopleWithSameBDay=3){

  bday_bins = seq( from=(1-0.5), to=(365+0.5), by=1 )
  nbins = length(bday_bins)-1
  
  bdMatch = matrix( 0, nrow=nSims, ncol=1 ) # 0/1 is there is a match or not
  for( ii in 1:nSims ){
    dbdays = sample( 1:365, nPeopleInRoom, replace=T ) # draw a room of people (and their bdays)
    hdb = hist( dbdays, breaks=bday_bins, plot=F )$counts
    if( any( hdb >= nPeopleWithSameBDay ) ){
      bdMatch[ii] = 1
    }
  }
   sum(bdMatch)/nSims
}

for( n in c(10,50,100,500,1000,5000,10000) ){
  print( sprintf( "nSims= %5d fraction with three or more of the same BDs=%10.6f", n, birthday_problem(nSims=n) ) )
}