Because Tymen isn’t familiar with the game craps…

• From www.online-craps-gambling.com/craps-rules.html

For the dice thrower (shooter) the object of the game is to throw a 7 or an 11 on the first roll (a win) and avoid throwing a 2, 3 or 12 (a loss). If none of these numbers (2, 3, 7, 11 or 12) is thrown on the first throw (the Come-out roll) then a Point is established (the point is the number rolled) against which the shooter plays. The shooter continues to throw until one of two numbers is thrown, the Point number or a Seven. If the shooter rolls the Point before rolling a Seven he/she wins, however if the shooter throws a Seven before rolling the Point he/she loses.

• So… describing set of permutations/possible outcomes of rolling two dice: (6,6) (6,5) (6,4) (6,3) (6,2) (6,1) (5,6) (5,5) (5,4) (5,3) (5,2) (5,1) (4,6) (4,5) (4,4) (4,3) (4,2) (4,1) (3,6) (3,5) (3,4) (3,3) (3,2) (3,1) (2,6) (2,5) (2,4) (2,3) (2,2) (2,1) (1,6) (1,5) (1,4) (1,3) (1,2) (1,1)

I promise I’ll make this much simpler

• Promised to get you something by morning
• for simulations 0=lose, 1=win
• some of this is patched together
• I promise I’ll remove the extraneous parts tomorrow after hopefully sleeping.

set.seed(888)#I'm just doing this so replications will produce same resuls
#just settting up a vector to be filled as we go
wins <- rep(0, 1000)
for(i in 1:1000){
#  simulate a Come-out rollif shooter wins on Come-out, 
wins[i] <- 1
}

makepoint <- function(point){
  made <- NA
  while(is.na(made)){ 
    roll <- sum(sample(6, 2, replace=T))
    if(roll == point){
      made <- TRUE
    } else if(roll == 7){
      made <- FALSE
    }
  } 
  return(made)
}
sum(wins) 

## [1] 1000

Tymen’s understanding of correctly simulating craps

simulating_craps <- function () {
  roll <- sum(sample(6,2,replace=T))
  if(roll==7 || roll==11)
    win <- T
  else if(roll==2 || roll==3 || roll==12)
    win <- F else
      win <- makepoint(roll)
    return(win)
}

• I looked up the Expected Value for Craps… E(X) =.493
• so I think my simulations are operating pretty accurately

n_sim <- 1000
wins <- 0
for(i in 1:n_sim)
wins <- wins + simulating_craps()
print(wins/n_sim) 

## [1] 0.486

finally what you actually want… here’s looking at likelihood

• note that I’m running it with with 4 separate simulation n’s
• but that’s just for curiousity sake.
• The resulting plot only uses the n=81 simulation

n_sim <- c(3, 9, 27, 81);
th <- seq(0, 1, length=200);
lik <- matrix(NA, 200, length(n_sim));

for (i in seq(along=n_sim)) 
{
  wins <- 0
   for(j in 1:n_sim[i])
   wins <- wins + simulating_craps()
   lik[,i] <- dbinom(wins, n_sim[i],th)
   lik[,i] <- lik[,i]/max(lik[,i])
}

Likelihood function for the parameter θ of winning a game of craps.

• emphasis on the word function:
• The likelihoods can be retrieved through calling the fourth simulation - in the matrix ‘lik’. ie: lik[,4]
• but probably want to decrease the number from 1000
• line 24 in the code

matplot(th, lik[,4], type="l", col=1, lty=1, 
xlab=expression(theta), ylab="likelihood", main="Likelihoods - 81 simulations")