3.5.1

Extend the simulation of throwing two dice in each of the following ways (one at a time, cumulatively):

a) Instead of 50 throws, extend the spreadsheet to make 500 row

Excel

R

Roll <- sample(1:6, 500, replace = T) + sample(1:6, 500, replace = T)

hist(Roll)

knitr::kable(describe(Roll))
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 500 6.946 2.457471 7 6.93 2.9652 2 12 10 0.0451817 -0.7416498 0.1099014

b)

Load the die by changing the probabilities of the faces to something other than uniform at 1/6 for each face. Be careful to ensure that your probabilties sum to 1.

Excel

R

Roll1 <- sample(1:6, 500, replace = T, prob = c(1/4, 1/4, 1/8, 1/8, 1/8, 1/8))
Roll2 <- sample(1:6, 500, replace = T, prob = c(1/4, 1/4, 1/8, 1/8, 1/8, 1/8))
Roll <- Roll1 + Roll2

hist(Roll)

knitr::kable(describe(Roll))
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 500 5.982 2.483145 6 5.8325 2.9652 2 12 10 0.4037795 -0.5346656 0.1110496

Simio

Modified based on the Monte Carlo Simulation solution provide in class:

c)

Use @RISK or another Excel add-in for static Monte Carlo spreadsheet simulation, to make 10,000 throws of the pair of fair dice, and compare your results to the true probabilities of getting the sum equal o each of 2, 3,….., 12 as well as the true expected value of 7.

Roll <- sample(1:6, 10000, replace = T) + sample(1:6, 10000, replace = T)
 
Probability <- c(1/36, 2/36, 3/36, 4/36, 5/36, 6/36, 5/36, 4/36, 3/36, 2/36, 1/36)
Simulated  <- table(Roll)/length(Roll)
Simulated Results
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 10000 6.9592 2.399856 7 6.957875 2.9652 2 12 10 0.0133225 -0.6088957 0.0239986
Probability Simulated
2 0.0277778 0.0292
3 0.0555556 0.0549
4 0.0833333 0.0806
5 0.1111111 0.1162
6 0.1388889 0.1464
7 0.1666667 0.1641
8 0.1388889 0.1377
9 0.1111111 0.1112
10 0.0833333 0.0796
11 0.0555556 0.0540
12 0.0277778 0.0261

3.5.5

In the Monte Carlo integration of Section 3.2.2, add to the spreadsheet calculation of the standard deviation of the 50 individual values, and use that, together with the mean already in cell H4, to compute a 95% confidence interval on the exact integral in cell I4.

Problem 3.5.17

Walther has a roadside produce stand where he sells oats, peas, beans, and barley. He buys these products at per-pound wholesale prices of, respectively, $1.05, $3.17, $1.99, and $0.95; he sells them at per-pound retail prices of, respectively, $1.29, $3.76, $2.33, and $1.65. Each day the amount demanded (in pounds) could be as little as zero for each product, and as much as 10, 8, 14, and 11 for oats, peas, beans, and barley, respectively; he sells only whole-pound amounts, no partial pounds. Assume a discrete uniform distribution for daily demand for each product over its range; assume as well that Walther always has enough inventory to satisfy all demand. The summer selling season is 90 days, and demand each day is independent of demand on other days. Create a spreadsheet simulation that will, for each day as well as for the whole season, simulate Walther’s total cost, total revenue and total profit.

crop_data <- 
  data.frame(
     crop = c("oats", "peas", "beans", "barley"),
     wac = c(1.05, 3.17, 1.99, 0.95),     # acquisition cost
     sell = c(1.29, 3.76, 2.23, 1.65),    # selling price
     max_demand = c(10, 8, 14, 11),       # maximum demand for this crop
     sim_demand = 0,                      # simulated daily demand for the given crop; uniform distribution
     cost = 0,                            # crop simulated cost (wac * simulated demand) 
     revenue = 0,                         # crop simulated revenue value (selling price * simulated demand)  
     profit = 0,                          # crop simulated profit value (selling price * simulated demand)
     stringsAsFactors = F) 
                        
         

sim_crop_season <- function(run_num, days) {
  
  # replicate each crop entry by the number of days in the season to simulate
  crop__sim <- as.data.frame(lapply(crop_data, rep, days))
  
  crop__sim$run_num <- run_num
  
  # set the simulated daily demand for the crop
  # This is a discrete uniform distribution for daily demand for each product over its range
  crop__sim$sim_demand <- round(runif(days, 0, crop__sim$max_demand), digits=0)  

  # calcualate the daily simulated cost, revenue, and profit for each crop  
  crop__sim <- 
               within(crop__sim,
                       { cost = wac * sim_demand
                         revenue = sell * sim_demand
                         profit = (sell * sim_demand) - (wac * sim_demand)
                       }
               )
  
}

Run the 90-day Simulation 100 times

sim <-  rbindlist(lapply(1:100, function(x) sim_crop_season(x, 90))) %>% as.data.frame()

Simulation results of a 90-day crop season (daily):

Simulated Crop Season Aggregate Metrics

sim %>% 
       group_by(run_num) %>% 
       summarise(cost    = sum(cost),
                 revenue = sum(revenue),
                  profit = sum(profit)) -> agg
  gather(agg, variable, value, -run_num) %>% arrange(run_num, variable) %>%
  ggplot( aes(x=variable, y = value, fill = variable)) +  geom_boxplot( alpha=0.7 ) +
          ggtitle("Crop Simulation Results") +
        theme_bw() + labs(x = "Metric", y="Total ($)") +
        theme(plot.title = element_text(size = 14, face = "bold"),
              text = element_text(size = 12),
              axis.title = element_text(face="bold"),
              axis.text.x=element_text(size = 11)) +
        scale_fill_brewer(palette = "Accent")
## Warning: package 'bindrcpp' was built under R version 3.3.3

Details