###########Example 1.2 from Lavine's book#################

#Set number of MC replications
N=1000

#make a vector of length N, filled with 0's
wins=rep(0, N)

#simulate a Come-out roll and if shooter wins on Come-out,
#set wins=1
for( i in 1:N) {
d <- sample( 1:6, 2, replace=T)
if(sum(d)==7 || sum(d)==11 ) wins[i]=1
}

#print the proportion of wins
sum(wins)/N 
#print true probability:
2/9

#Execute the above as an R function:

prob.win=function(N=1000){
wins=rep(0, N)
for(i in 1:N){
d <- sample( 1:6, 2, replace=T)
if(sum(d)==7 || sum(d)==11 ) wins[i]=1
}
cat(c("true prob=",round(2/9,4)),fill=T)
cat(c("empirical prob=",round(mean(wins),4),"[using", N, "replications]"),fill=T)
}

prob.win(N=100)
prob.win(N=100)
prob.win(N=1000)
prob.win(N=1000)
prob.win(N=10000)
prob.win(N=10000)
prob.win(N=100000)
prob.win(N=100000)


#Compute the probability of loss:
prob.loss=function(N=1000){
loss=rep(0,N)
for(i in 1:N){
d <- sample( 1:6, 2, replace=T)
if(sum(d)==2 || sum(d)==3 || sum(d)==12 ) loss[i]=1
}
cat(c("true prob=",round(1/9,4)),fill=T)
cat(c("empirical prob=",round(mean(loss),4),"[using", N, "replications]"),fill=T)
}

prob.loss(N=100)
prob.loss(N=100)
prob.loss(N=1000)
prob.loss(N=1000)
prob.loss(N=10000)
prob.loss(N=10000)
prob.loss(N=100000)
prob.loss(N=100000