Chapter 6 Homework

Kelton Chapter 6 problems

1. Use file from text to fit one or more probability distributions, including goodness-of-fit testing and probability plots. What’s the recommendation for a distribution to be used in the simulation model from which to generate these inter-arrival times? Provide the correct Simio expression, heading any parameterization-difference issues.

library(xlsx)
prob.06.01.df <- read.xlsx("~./GitHub/MSDA---Coursework/Data 604/Week 9 - HW/Problem_Dataset_06_01.xls", sheetIndex = 1)

library(fitdistrplus)
plotdist(prob.06.01.df$X27.88, histo = TRUE)

descdist(prob.06.01.df$X27.88, boot = 100)

summary statistics
------
min:  2.06   max:  48.65 
median:  8.74 
mean:  11.53122 
estimated sd:  9.621378 
estimated skewness:  1.968746 
estimated kurtosis:  7.88031 

Based on our Cullen and Frey Graph the exponential distribution is our best fit.

dis <- fitdist(prob.06.01.df$X27.88, distr = "exp")
summary(dis)

Fitting of the distribution ' exp ' by maximum likelihood 
Parameters : 
      estimate Std. Error
rate 0.0867211 0.01354176
Loglikelihood:  -141.2474   AIC:  284.4948   BIC:  286.2083 

plot(dis)

In Simio it would be “Random.Exponential(0.0867211)”

2. Use file from text to fit one or more probability distributions, including goodness-of-fit testing and probability plots. What’s the recommendation for a distribution to be used in the simulation model from which to generate these call-duration times? Provide the correct Simio expression, heading any parameterization-difference issues.

library(fitdistrplus)
prob.06.02.df <- read.xlsx("~./GitHub/MSDA---Coursework/Data 604/Week 9 - HW/Problem_Dataset_06_02.xls", sheetIndex = 1)
plotdist(prob.06.02.df$X36.13, histo = TRUE)

descdist(prob.06.02.df$X36.13, boot = 100)

summary statistics
------
min:  3.8   max:  27.86 
median:  8.91 
mean:  10.3437 
estimated sd:  5.191441 
estimated skewness:  1.390486 
estimated kurtosis:  5.02623 

dis <- fitdist(prob.06.02.df$X36.13,  "gamma")
dis

Fitting of the distribution ' gamma ' by maximum likelihood 
Parameters:
       estimate Std. Error
shape 4.8173431  0.9716745
rate  0.4657621  0.0990183

plot(dis)

coef(dis)

    shape      rate 
4.8173431 0.4657621 

In Simio its Random.Gamma( 4.8173431 , 0.4657621 )

3. Use file from text to fit one or more probability distributions, including goodness-of-fit testing and probability plots. What’s the recommendation for a distribution to be used in the simulation model from which to generate the number of extra tech-support people need for a call? Provide the correct Simio expression, heading any characterization-difference issues.

prob.06.03.df <- read.xlsx("~./GitHub/MSDA---Coursework/Data 604/Week 9 - HW/Problem_Dataset_06_03.xls", sheetIndex = 1)

library(fitdistrplus)
plotdist(prob.06.03.df$X2, histo = TRUE)

descdist(prob.06.03.df$X2, boot = 100)

summary statistics
------
min:  1   max:  8 
median:  3 
mean:  2.977273 
estimated sd:  1.547521 
estimated skewness:  0.9862332 
estimated kurtosis:  4.349472 

dis <- fitdist(prob.06.03.df$X2, distr = "lnorm")
summary(dis)

Fitting of the distribution ' lnorm ' by maximum likelihood 
Parameters : 
         estimate Std. Error
meanlog 0.9541848 0.08153836
sdlog   0.5408643 0.05765544
Loglikelihood:  -77.37561   AIC:  158.7512   BIC:  162.3196 
Correlation matrix:
        meanlog sdlog
meanlog       1     0
sdlog         0     1

plot(dis)

coef(dis)

  meanlog     sdlog 
0.9541848 0.5408643 

In Simio its Random.Lognormal(0.9541848, 0.5408643)

4. Derive the inverse-CDF formula for generation random variants from a continuous uniform distribution between real numbers a and b (a < b )

CDF

\[f(x) = \int_a^{b} \frac{x-a}{b-a} dx\] Inverse

\[f^{-1}(x) = a + (b - a) x\]

5. Derive the inverse-CDF formula for generation random variate from a Weibull distribution. Look up its definition (including its CDF) in one of the references, either book or online, cited in this chapter. Check the simio definition in its documentation to be sure you have your formula parameterized proper.

\[f(x) = 1 - e^{1-(1/\alpha)x^\beta} \] \[f^{-1}(x)= (-\alpha ln(1-u))^{1/\beta}\]

8. Problem in text:

library(tidyverse)

Loading tidyverse: ggplot2
Loading tidyverse: tibble
Loading tidyverse: tidyr
Loading tidyverse: readr
Loading tidyverse: purrr
Loading tidyverse: dplyr
package <U+393C><U+3E31>tidyr<U+393C><U+3E32> was built under R version 3.4.2package <U+393C><U+3E31>purrr<U+393C><U+3E32> was built under R version 3.4.2package <U+393C><U+3E31>dplyr<U+393C><U+3E32> was built under R version 3.4.2Conflicts with tidy packages -----------------------------------------------------------------------------------------------------------------
filter(): dplyr, stats
lag():    dplyr, stats
select(): dplyr, MASS

produce <- NULL
produce$oats <-  list(cost = 1.05, 
                      revenue =  1.29,
                      max_demand = 10, 
                      sizes = c(0, .5 , 1, 1.5, 2, 3, 4, 5, 7.5, 10),
                      prob = c(.05, .07, .09, .11, .15, .25, .1, .09, .06, .03)
                      )
produce$peas <-  list(cost = 3.17, 
                      revenue =  3.76,
                      max_demand = 8, 
                      sizes = c(0, .5, 1, 1.5, 2, 3),
                      prob = c(.1, .2, .2, .3, .1, .1)
                      )
produce$beans <-  list(cost = 1.99, 
                      revenue =  2.23,
                      max_demand = 14, 
                      sizes = c(0, 1, 3, 4.5),
                      prob = c(.2, .4, .3, .1)
                      )
produce$barley <- list(cost = .95, 
                      revenue =  1.65,
                      max_demand = 11, 
                      sizes = c(0, .5, 1, 3.5),
                      prob = c(.2, .4, .3, .1)
                      )
for (i in names(produce)){
  produce[[i]]["cumprob"] <- list(cumsum(produce[[i]]$prob))
}

library(data.table)
for (i in names(produce)){
  x <- NULL
  for (j in seq(1,90)){
  x <- append(x, round(runif(1),4))}
  produce[[i]]["observedprob"] <- list(x)
  x <- NULL
  for (t in seq(1,90)){
  x <- append(x,(produce[[i]]$sizes[sum((produce[[i]]$cumprob - produce[[i]]$prob) <= produce[[i]]$observedprob[[t]])]))
  }
  produce[[i]]["lbs.sold"]   <- list(x)
  produce[[i]]["total.cost"] <- list(produce[[i]]$lbs.sold * produce[[i]]$cost)
  produce[[i]]["total.rev"]  <- list(produce[[i]]$lbs.sold * produce[[i]]$revenue)
}
for (i in names(produce)){
print(paste(i, "total cost = ", sum(produce[[i]]$total.cost)))
print(paste(i, "total revenue = ", sum(produce[[i]]$total.rev)))
print(paste(i, "profit = ", sum(produce[[i]]$total.rev - produce[[i]]$total.cost)))
  
}

[1] "oats total cost =  265.65"
[1] "oats total revenue =  326.37"
[1] "oats profit =  60.72"
[1] "peas total cost =  342.36"
[1] "peas total revenue =  406.08"
[1] "peas profit =  63.72"
[1] "beans total cost =  320.39"
[1] "beans total revenue =  359.03"
[1] "beans profit =  38.64"
[1] "barley total cost =  68.4"
[1] "barley total revenue =  118.8"
[1] "barley profit =  50.4"