#
# 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 31:
#
mixed_up_mail = function(nSims=10,nEnvelopes=10){

  bins = seq( from=(0-0.5), to=(nEnvelopes+0.5), by=1 )

  nMatches = matrix( 0, nrow=nSims, ncol=1 )
  for( ii in 1:nSims ){
    samp = sample( 1:nEnvelopes, nEnvelopes )
    nMatches[ii] = sum( samp == 1:nEnvelopes )
  }

  hist( nMatches, breaks=bins, plot=FALSE )$intensities
}

nEnvelopes = 10

nDraws = c(10,100,1000)

M = matrix( 0, nrow=nEnvelopes+1, ncol=length(nDraws) )
for( ii in 1:length(nDraws) ){
  mum = mixed_up_mail( nSims=nDraws[ii] )
  M[,ii] = mum
}

#postscript("../../WriteUp/Graphics/Chapter2/prob_31_hist.eps", onefile=FALSE, horizontal=FALSE)

barplot( t(M), beside=T, main="MC distribution", ylab="frequency", xlab="10,100,1000 number of experiments" )

#dev.off()

# whats the expected number of rolls needed to get three looks like it is ~ 2
#
t(M) %*% matrix( 1:(length(bins)-1), nrow=(length(bins)-1), ncol=1 )