Homework 5 Probability Models

Problem 1 — Estimating Auto Accident Costs

# Problem1_R_simulation.R
set.seed(12345)
n.sim <- 5000
n.trips <- 500
p_acc <- 0.001
alpha <- 1.5
beta <- 3
min_cost <- 500
max_cost <- 20000

total_costs <- numeric(n.sim)

for (i in 1:n.sim) {
  k <- rbinom(1, size = n.trips, prob = p_acc)
  if (k == 0) {
    total_costs[i] <- 0
  } else {
    # sample k Beta(1.5,3) then scale
    x <- rbeta(k, shape1 = alpha, shape2 = beta)
    costs <- min_cost + x * (max_cost - min_cost)
    total_costs[i] <- sum(costs)
  }
}

# Summary stats
mean_total <- mean(total_costs)
sd_total <- sd(total_costs)
# 95% CI for mean (normal approximation using sample mean and sd/sqrt(n.sim))
se_mean <- sd_total / sqrt(n.sim)
ci_mean <- mean_total + c(-1,1) * qnorm(0.975) * se_mean

# Probability total cost > 8000 with 95% CI (binomial on sims)
prop_gt8000 <- mean(total_costs > 8000)
# CI for proportion using Wilson (prop.test gives approx)
prop_test <- prop.test(sum(total_costs > 8000), n.sim)
ci_prop <- prop_test$conf.int

# Print results
cat("Mean total cost:", round(mean_total,2), "\n")

## Mean total cost: 3466.36

cat("SD total cost:", round(sd_total,2), "\n")

## SD total cost: 5643.95

cat("95% CI for mean total cost:", round(ci_mean,2), "\n")

## 95% CI for mean total cost: 3309.92 3622.8

cat("P(total cost > 8000):", round(prop_gt8000,4), "\n")

## P(total cost > 8000): 0.1918

cat("95% CI for that proportion:", round(ci_prop,4), "\n")

## 95% CI for that proportion: 0.181 0.2031

Problem 2 — Project Management

# Problem2_R_simulation.R
set.seed(23456)
n.sim <- 5000

# Define the discrete distributions as vectors and probabilities
A_vals <- c(5,6,7,8); A_p <- c(0.25, 0.35, 0.25, 0.15)
B_vals <- c(3,5,7);       B_p <- c(0.20, 0.55, 0.25)
C_vals <- c(10,12,14,16,18); C_p <- c(0.10,0.25,0.40,0.20,0.05)
D_vals <- c(8,10);        D_p <- c(0.60, 0.40)

proj_length <- numeric(n.sim)

for (i in 1:n.sim) {
  a <- sample(A_vals, size=1, prob=A_p)
  b <- sample(B_vals, size=1, prob=B_p)
  c_ <- sample(C_vals, size=1, prob=C_p)
  d <- sample(D_vals, size=1, prob=D_p)
  proj_length[i] <- a + b + c_ + d
}

mean_proj <- mean(proj_length)
sd_proj <- sd(proj_length)
# probability completed in 35 weeks or less
prob_le_35 <- mean(proj_length <= 35)
# 95% CI for that probability using prop.test
pt <- prop.test(sum(proj_length <= 35), n.sim)
ci_prob <- pt$conf.int

cat("Mean project length:", round(mean_proj,3), "\n")

## Mean project length: 33.91

cat("SD project length:", round(sd_proj,3), "\n")

## SD project length: 2.82

cat("P(<=35 weeks):", round(prob_le_35,4), "\n")

## P(<=35 weeks): 0.7206

cat("95% CI for P(<=35):", round(ci_prob,4), "\n")

## 95% CI for P(<=35): 0.7079 0.733

Problem 3 — Wearable Product Launch

# Problem3_R_simulation.R
set.seed(34567)
n.sim <- 5000
fixed_cost <- 350000
price <- 200
dmean <- 9000
dsd <- 2000

profits <- numeric(n.sim)
units_sold_vec <- integer(n.sim)
variable_costs <- numeric(n.sim)

for (i in 1:n.sim) {
  # variable cost per unit
  var_cost <- runif(1, min = 90, max = 150)
  # demand (round and truncate at 0)
  demand <- round(rnorm(1, mean = dmean, sd = dsd))
  if (demand < 0) demand <- 0
  revenue <- price * demand
  total_var_cost <- var_cost * demand
  profit <- revenue - fixed_cost - total_var_cost
  profits[i] <- profit
  units_sold_vec[i] <- demand
  variable_costs[i] <- var_cost
}

mean_profit <- mean(profits)
sd_profit <- sd(profits)
prob_loss <- mean(profits < 0)
pt3 <- prop.test(sum(profits < 0), n.sim)
ci_loss <- pt3$conf.int

cat("Mean profit:", round(mean_profit,2), "\n")

## Mean profit: 369150.5

cat("SD profit:", round(sd_profit,2), "\n")

## SD profit: 228850.1

cat("P(loss):", round(prob_loss,4), "\n")

## P(loss): 0.0302

cat("95% CI for P(loss):", round(ci_loss,4), "\n")

## 95% CI for P(loss): 0.0257 0.0354