> #### Simulate 20 rolls of a die using 'sample' function
>
> sim1.out<-sample(c(1,2,3,4,5,6), 20, replace=T)
> sim1.out
 [1] 1 1 6 5 2 5 2 2 5 3 5 2 5 6 2 1 2 2 4 1
> 
> # Specify probabilities
> sim2.out<-sample(c(1,2,3,4,5,6), 20, replace=T, prob=c(1/9, 1/9, 1/9, 1/9, 2/9, 3/9))
> sim2.out
 [1] 2 5 6 6 3 6 6 4 6 5 4 1 5 2 6 3 6 2 6 6
> 
> # Calculate the ratio Y/X where Y is bigger and X is smaller of the two rolls
> roll1 <- sample(c(1,2,3,4,5,6), 1)
> roll2 <- sample(c(1,2,3,4,5,6), 1)
> Y<-max(roll1,roll2)
> X<-min(roll1,roll2)
> Z<-Y/X
> cbind(roll1, roll2, Y, X, Z)
     roll1 roll2 Y X Z
[1,]      4      1 4 1 4
> 
> ######## SIMULATIONS: 'REPLICATE' set of commands n=50 times, store output
>
> store1<-replicate(50, {
+ roll1 <- sample(c(1,2,3,4,5,6), 1)
+ roll2 <- sample(c(1,2,3,4,5,6), 1)
+ Y<-max(roll1,roll2)
+ X<-min(roll1,roll2)
+ Z<-Y/X
+ })
> store1
 [1] 1.500000 1.666667 1.500000 1.200000 4.000000 1.250000 1.500000 1.250000 2.000000 1.666667 3.000000 2.000000 1.666667
[14] 1.500000 2.000000 4.000000 1.000000 2.500000 2.000000 2.000000 2.500000 6.000000 2.000000 2.500000 6.000000 1.000000
[27] 5.000000 5.000000 1.000000 3.000000 2.000000 1.000000 1.200000 2.000000 1.666667 2.000000 2.000000 2.000000 1.666667
[40] 1.333333 1.200000 6.000000 1.333333 1.500000 2.000000 3.000000 1.500000 2.000000 2.000000 5.000000

> ######## SIMULATIONS USING 'FOR' LOOP
>
> store2 <-rep(NA, 50)   # Initialize storage 
> for (i in 1:50) {
+ roll1 <- sample(c(1,2,3,4,5,6), 1)
+ roll2 <- sample(c(1,2,3,4,5,6), 1)
+ store2[i]<-max(roll1,roll2)/min(roll1,roll2)
+ }
> store2
 [1] 3.000000 1.333333 2.000000 1.500000 2.500000 2.000000 4.000000 3.000000 1.000000 4.000000 2.500000 3.000000 1.200000
[14] 1.666667 3.000000 2.500000 1.000000 3.000000 1.666667 1.333333 1.500000 1.000000 1.250000 5.000000 1.000000 1.333333
[27] 1.500000 1.000000 1.500000 1.200000 1.500000 1.333333 2.000000 6.000000 1.333333 1.333333 1.333333 3.000000 1.666667
[40] 2.000000 5.000000 1.500000 1.250000 2.000000 1.000000 1.666667 1.333333 2.500000 1.500000 1.000000

> ######## SIMULATIONS USING VECTORIZATION
>
> vec1<-sample(c(1,2,3,4,5,6), 50, replace=T)
> vec2<-sample(c(1,2,3,4,5,6), 50, replace=T)
> vec1
 [1] 6 5 5 6 4 2 1 1 1 2 6 4 6 6 4 3 2 4 6 6 3 1 2 6 1 2 3 1 3 2 6 5 2 3 1 2 5 1 5 2 3 6 1 2 6 4 3 3 4 3
> vec2
 [1] 3 4 5 3 1 5 2 1 1 1 6 1 6 1 3 1 3 2 5 3 3 6 5 5 1 6 4 4 2 1 2 6 1 1 6 3 2 3 4 2 5 5 1 1 5 4 2 4 4 3
> store3<-apply(cbind(vec1,vec2),1,max)/apply(cbind(vec1,vec2),1,min)
> store3
 [1] 2.000000 1.250000 1.000000 2.000000 4.000000 2.500000 2.000000 1.000000 1.000000 2.000000 1.000000 4.000000 1.000000 6.000000 1.333333 3.000000
[17] 1.500000 2.000000 1.200000 2.000000 1.000000 6.000000 2.500000 1.200000 1.000000 3.000000 1.333333 4.000000 1.500000 2.000000 3.000000 1.200000
[33] 2.000000 3.000000 6.000000 1.500000 2.500000 3.000000 1.250000 1.000000 1.666667 1.200000 1.000000 2.000000 1.200000 1.000000 1.500000 1.333333
[49] 1.000000 1.000000