Newton-Pepys ver2

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

Run Sequentially

# Stat 6401
# Example 1.4.19
# Newton-Pepys problem
B <- 1000000  # number of replications of the experiment
# A : At least one 6 appears when 6 fair dice are rolled.
n <- 6  # number of dice
k <- 1  # number of time 6 appears in the n dice
r <- replicate(B,sum(sample(1:6, n, replace = TRUE) == 6))
sum(r >= k)/B  # estimate probability

[1] 0.66503

# B : At least two 6's appears when 12 fair dice are rolled.
n <- 12  # number of dice
k <- 2  # number of time 6 appears in the n dice
r <- replicate(B,sum(sample(1:6, n, replace = TRUE) == 6))
sum(r >= k)/B  # estimate probability

[1] 0.618667

# C : At least two 6's appears when 18 fair dice are rolled.
n <- 18  # number of dice
k <- 3  # number of time 6 appears in the n dice
r <- replicate(B,sum(sample(1:6, n, replace = TRUE) == 6))
sum(r >= k)/B  # estimate probability

[1] 0.596521

As the number of dice increases the chance of winning decreases.

Run using a sequential for loop

# using loop
#start time
strt <- Sys.time()
games <- 3  # number of times 6 appears
B <- 1000000  # number of replications of the experiment
n <- 6  # number of dice
r <- numeric(games)
for(k in 1:games){
  r[k] <- sum(replicate(B,sum(sample(1:6, n*k, replace = TRUE) == 6)) >= k)/B
}
prob.est <- r
prob.est

[1] 0.664172 0.619559 0.597354

#end time
print(Sys.time()-strt)

Time difference of 27.44176 secs

Run in parallel

# parallel code
#start time
strt <- Sys.time()
games <- 3  # number of times 6 appears
library(foreach)
library(doParallel)
cl <- makeCluster(games)
registerDoParallel(cl)
B <- 1000000  # number of replications of the experiment
n <- 6  # number of dice
r <- numeric(games)
prob.est <- foreach(k = 1:games, .combine=rbind) %dopar% {
  r[k] <- sum(replicate(B,sum(sample(1:6, n*k, replace = TRUE) == 6)) >= k)/B
}
prob.est

             [,1]
result.1 0.665022
result.2 0.617407
result.3 0.597322

stopCluster(cl)
#end time
print(Sys.time()-strt)

Time difference of 12.46573 secs

# parallel code for linux or OSX
#start time
strt <- Sys.time()
games <- 3  # number of times 6 appears
library(foreach)
library(doMC)
registerDoMC(games) 
B <- 1000000  # number of replications of the experiment
n <- 6  # number of dice
r <- numeric(games)
prob.est <- foreach(k = 1:games, .combine=rbind) %dopar% {
  r[k] <- sum(replicate(B,sum(sample(1:6, n*k, replace = TRUE) == 6)) >= k)/B
}
prob.est

             [,1]
result.1 0.664976
result.2 0.618620
result.3 0.597547

#end time
print(Sys.time()-strt)

Time difference of 10.65973 secs

# parallel code for Windows
#start time
strt <- Sys.time()
games <- 3  # number of times 6 appears
library(foreach)
library(doSNOW)
cl <- makeCluster(games)
B <- 1000000  # number of replications of the experiment
n <- 6  # number of dice
r <- numeric(games)
prob.est <- foreach(k = 1:games, .combine=rbind) %dopar% {
  r[k] <- sum(replicate(B,sum(sample(1:6, n*k, replace = TRUE) == 6)) >= k)/B
}
prob.est

             [,1]
result.1 0.665028
result.2 0.618497
result.3 0.597397

stopCluster(cl)
#end time
print(Sys.time()-strt)

Time difference of 11.41324 secs

Using parallel processing reduces the overall time to run the simulations.

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.