SimulationExam1

Problem 1: Yahtzee

Part A.

Write a function that simulates the rolling of 5 six-sided die. This function should return a vector of length 5.

dice <- function(n){
sample(1:6, size = 5, replace = TRUE)
}
x <- dice()
x

## [1] 4 3 2 1 4

Part B.

A small straight occurs when four of the dice are sequential numbers. Using the function written in the previous part, estimate the probability using Monte Carlo simulation that you will roll a small straight on a SINGLE roll of the five dice (i.e. You roll all five dice once, and exactly four of them are sequential).

A small straight is four consecutive numbers, not necessarily in order. The three possible small straights are 1, 2, 3, 4, 2, 3, 4, 5, and 3, 4, 5, 6.

straight <- function(){
x <- dice()
if (all(c(1,2,3,4) %in% x) |
   all(c(2,3,4,5) %in% x) | 
  all(c(3,4,5,6) %in% x)) {
return(TRUE)
}  
else {
  return(FALSE)
}
}

straight()

## [1] FALSE

straightballin <- replicate(10000, straight())
mean(straightballin)

## [1] 0.1569

#Check
960/7776

## [1] 0.1234568

Part C

full <- function(){
  x <- dice()
  x <- table(x)
if (all(c(2,3) %in% x) | all(c(3,2) %in% x)){
  return(TRUE)
}
  else {
    return(FALSE)
}
}

full()  

## [1] FALSE

fullballin <- replicate(10000, full())
mean(fullballin)

## [1] 0.0401

#Check
(((5*4)/2)*30)/7776

## [1] 0.03858025

Problem D

nomatch <- function(){
  x <- dice()
  x <- table(x)
if ( (all(c(2) %in% x)) | (all(c(3) %in% x)) | (all(c(4) %in% x)) | (all(c(5) %in% x)) ){
  return(FALSE)
}
  else {
    return(TRUE)
}
}

nomatch()  

## [1] FALSE

nomatchbud <- replicate(10000, nomatch())
mean(nomatchbud)

## [1] 0.0955

#Check
(6*5*4*3*2)/7776

## [1] 0.09259259

Problem 2: Secant Method

I used code from the last homework.

From the plot we see there are many roots for the function. I focused on the negative roots and checked them in wolfram alpha as real roots.

func <- function(x){
  sin(exp(x)) + cos(exp(-x))
}

curve(expr = func, col='purple', xlim=c(-4,5))
abline(h=0)
abline(v=0)

secant.method <- function(f, x0, x1, tol = 1e-9, n = 500) {
  for (i in 1:n) {
    x2 <- x1 - f(x1) / ((f(x1) - f(x0)) / (x1 - x0)) 
    if (abs(x2 - x1) < tol) { 
      return(x2)
    }
    x0 <- x1
    x1 <- x2 
  }
}

secant.method(func, -4, -3)

## [1] -3.286195

#-3.286

secant.method(func, -3.5, -3)

## [1] -3.394897

#-3.394

secant.method(func, -3, -2.5)

## [1] -3.157827

#-3.157

secant.method(func, -2, -1)

## [1] -1.501776

#-1.501776

secant.method(func, -1, 0)

## [1] -0.7212255

#-0.7212255

Problem 3: Mersenne Twister Algorithm

This code was found in your notes. We mod the output by

merstwist <- function(a=1117,b=133,m=10){
  seed <- runif(1,0,1)
  out <- (a*seed + b) %%m
  return(out%%1)
}

merstwist(1117,133,10)

## [1] 0.1099677

Problem 4: Inversion Method

First invert the function and then use the above function.

frechet <- function(a){
  return( (-log(merstwist())) ^ (-1/a) )
}

frechet(1)

## [1] 0.8738638

Problem 5

Based on the previous problem, fix α = 1, sample from the distribution 1000 times and store this data. Repeat this with values of α = 2, 3, 4, 5, and 6. Create a 2 x 3 grid displaying the results of these random draws (i.e. the upper left plot should be the histogram corresponding to α = 1, upper right should correspond to α = 2, etc.) Make sure to label your grid.

par(mfrow=c(3,2))

n <- 1000
z <- rep(NA,n)
for (i in 1:n){
  z[i] <- frechet(1)
}
hist(z, xlim = c(0,6), main = paste("alpha = 1"))

n <- 1000
z <- rep(NA,n)
for (i in 1:n){
  z[i] <- frechet(2)
}
hist(z, xlim=c(0,6), main = paste("alpha = 2"))

n <- 1000
z <- rep(NA,n)
for (i in 1:n){
  z[i] <- frechet(3)
}
hist(z, xlim=c(0,6), main = paste("alpha = 3"))

n <- 1000
z <- rep(NA,n)
for (i in 1:n){
  z[i] <- frechet(4)
}
hist(z, xlim=c(0,6), main = paste("alpha = 4"))

n <- 1000
z <- rep(NA,n)
for (i in 1:n){
  z[i] <- frechet(5)
}
hist(z, xlim=c(0,6), main = paste("alpha = 5"))

n <- 1000
z <- rep(NA,n)
for (i in 1:n){
  z[i] <- frechet(6)
}
hist(z, xlim=c(0,6), main = paste("alpha = 6"))