Q1

Suppose a fraction 17% of the microchips produced by a leading manufacturer is defective. For the manufacturer, the number of microchips delivered to its client depends on how many microchips are inspected and accepted by 50 inspectors in the factory. Historically, given that a microchip is defective, an inspector wrongly accepts the chip 10% of the time, thinking it has no defect. If a microchip is not defective, an inspector, however, wrongly rejects the chip 5% of the time. In general, each of the 50 inspectors can inspect 20 chips per hour in a working day of 8 hours (lunch break NOT included in the 8 hours). Use the rbinom( ) and whatever functions needed in R to write a simulation program for the case above. In a working day, the microchips delivered to the client contain some good ones (correctly accepted) and some bad ones. In your simulation program, calculate the ratio of good ones to the sum of good and bad ones. The simulation program should be a function called dailychips. The function must return 1) the daily number of delivered chips (that is random) and 2) the daily ratio of good ones mentioned above. Simulate the daily operations for 1,000 runs and answer the following questions based on simulation results.

defect=0.17
notDefect=1-defect

defect.accepted=0.1
defect.notAccepted=1-defect.accepted

notDefect.notAccepted=0.05
notDefect.accepted=1-notDefect.notAccepted

num.inspector = 50
num.inspection = 20
inspector.workhr = 8


dailychips=function(n=num.inspector*num.inspection*inspector.workhr){
  acceptornot=rep(NA, n)
  chip.status=sample(c(1, 0),  # 1: no defect, 0: defect
                   n, 
                   prob=c(notDefect, defect), 
                   replace=TRUE)
  
  if(any(chip.status==1)){
    acceptornot[which(chip.status==1)]=
      rbinom(n=sum(chip.status==1), size=1, prob=notDefect.accepted)  # 1: accept, 0: not accept
  }
  
  if(any(chip.status==0)){
    acceptornot[which(chip.status==0)]=
      rbinom(n=sum(chip.status==0), size=1, prob=defect.accepted)
  }
  #return(cat("daily number of delivered chips:", sum(acceptornot), "\n the daily ratio of good ones", sum(chip.status[which(acceptornot==1)]==1) / sum(acceptornot)))
  return(c(sum(acceptornot), sum(chip.status[which(acceptornot==1)]==1) / sum(acceptornot)))
}

S = 1000
dailychips.table=replicate(S, dailychips())
dim(dailychips.table)
## [1]    2 1000

Q1-a

  1. The manufacturer claims that the everyday chips delivered to its client has at least 98% good ones. What is the probability that such a claim is true?
sum(dailychips.table[2,]>=0.98)/S
## [1] 0.272

Q1-b

  1. What is the probability that the manufacturer can deliver 6,400 microchips accepted by its inspectors in a typical working day of 8 hours?
sum(dailychips.table[1,]>=6400)/length(dailychips.table[1,])  # 有 6400 以上就表示也可以只 deliver 6400 ?
## [1] 0.9

Q2

Use the dailychips function in Q1 to simulate montlychips and quarterlychips. The monthlychips is the sum of delivered chips over 30 days, whereas the quarterlychips is the sum of delivered chips over 90 days respectively. Ignore good ratios in the two cases and focus on the number of delivered chips over 30 and 90 days. Generate 1,000 random samples for the two numbers (montly and quarterly) of delivered chips.

Generate the histograms (use hist( ) in R and set breaks=20) of the simulated montlychips and quarterlychips. Do they look like normal distributions?

Use the shaprio.test ( ) in R to perform normality tests for montlychips and quarterlychips (if p-value<=0.05, reject H0: Normality holds).

Do you find any evidence for normality? If yes, can you provide theoretical explanations for that (hint: check chapter 3 for a theorem)?

monthlychips = function(days=30) {
  table = replicate(days, dailychips())
  return(sum(table[1,]))
}

monthlychips.table=replicate(S, monthlychips())

hist(monthlychips.table,
     main = "Distribution of number of chips delivered monthly",
     xlab = "Number of chips delivered monthly",
     breaks = 20)

quarterlychips.table=replicate(S, monthlychips(days=90))

hist(quarterlychips.table,
     main = "Distribution of number of chips delivered quarterly",
     xlab = "Number of chips delivered quarterly",
     breaks = 20)

Shaprio.test

shapiro.test(monthlychips.table)
## 
##  Shapiro-Wilk normality test
## 
## data:  monthlychips.table
## W = 0.99795, p-value = 0.2614
shapiro.test(quarterlychips.table)
## 
##  Shapiro-Wilk normality test
## 
## data:  quarterlychips.table
## W = 0.99858, p-value = 0.6093

p-value 未小於0.05，沒有拒絕虛無假設，兩者皆符合常態分佈

The Central Limit Theorem，中央極限定理，無論原本變數產生的分佈為何，只要次數夠大(n>30)，最終的分佈會接近常態分佈

Q3

This case is adapted from Chapter 5 in Bertsimas & Freund (2004).

profit.permonth = function(months=1){
  restaur.seat = 50
  cost.nonlaborfixed.permonth = 3995
  cost.eachmeal = 11
  cost.labor.permonth=runif(1, 5040, 6860)
  cost.permonth = cost.nonlaborfixed.permonth + cost.labor.permonth
  
  num.mealsold.permonth=rnorm(1, 3000, 1000)
  price.eachmeal=sample(c(20, 18.5, 16.5, 15), 1, c(0.25, 0.35, 0.3, 0.1), replace=TRUE)
  revenue.eachmeal = price.eachmeal - cost.eachmeal
  revenue.permonth = sum(revenue.eachmeal*num.mealsold.permonth)
  
  profit.permonth = revenue.permonth - cost.permonth
  return(profit.permonth)
}

salary.restaur.withpartner = function() {
  profit.permonth = profit.permonth()
  if (profit.permonth < 3500) {
    salary.restaur.withpartner = 3500
  }
  else if (profit.permonth >= 3500 && profit.permonth <= 9000) {
    salary.restaur.withpartner = profit.permonth
  }
  else {
    salary.restaur.withpartner = 9000 + (profit.permonth - 9000)*0.1
  }
  return(salary.restaur.withpartner)
}

Q3-1

Without considering the partnership opportunity, what would be Sanjay’s expected monthly salary at Gentle Lentil? How does this compare to his monthly salary at the consulting firm?

S = 1000
profit.permonth.table=replicate(S, profit.permonth())
cat("Sanjay’s expected monthly salary at Gentle Lentil (without partnership) will be ", sum(profit.permonth.table)/S)
## Sanjay’s expected monthly salary at Gentle Lentil (without partnership) will be  10697.34
salary.job.peryear = 80000
salary.job.permonth = salary.job.peryear/12
cat("Sanjay’s expected monthly salary at the consulting firm will be ", salary.job.permonth)
## Sanjay’s expected monthly salary at the consulting firm will be  6666.667

Q3-2

With considering the partnership opportunity, what would be Sanjay’s expected monthly income at Gentle Lentil? How does this compare to his monthly salary at the consulting firm and his salary without the partnership at Gentle Lentil?

S = 1000
salary.restaur.withpartner.table =replicate(S, salary.restaur.withpartner())
cat("Sanjay’s expected monthly salary at Gentle Lentil (with partnership) will be ", sum(salary.restaur.withpartner.table)/S)
## Sanjay’s expected monthly salary at Gentle Lentil (with partnership) will be  7707.43

Q3-3

If you were Sanjay, what quantitative questions about the potential salary at Gentle Lentil (with and without the partnership opportunity) would you want to answer with the aid of simulation modeling, before deciding whether to launch the Gentle Lentil Restaurant or to accept the consulting position? Show answers to the quantitative questions based on your simulation model.

##### 開餐廳 (with partnership) 的收入 高於 consulting firm salary 的機率：
sum(salary.restaur.withpartner.table>6666.7)/S
## [1] 0.667
##### 開餐廳 (without partnership) 的收入 高於 consulting firm salary 的機率：
sum(profit.permonth.table>6666.7)/S
## [1] 0.651