# Discrete Random Variables
# Objective: Map outcomes to one integer-valued random variable and build its support.

# Module 6 R Note 1: Random variables
outcomes <- c("Low demand", "Typical demand", "High demand")
values <- c(40, 55, 70)
mapping <- data.frame(outcome = outcomes, X = values)
stopifnot(length(unique(mapping$outcome)) == nrow(mapping))
print(mapping)
print(sort(unique(mapping$X)))
print("X is daily demand; x is one possible demand value.")

# ----

# Distribution Validation
# Objective: Validate a custom PMF without silently repairing it.

# Module 6 R Note 2: Validate a PMF
demand <- read.csv("module6_custom_discrete_demand_synthetic.csv")
pmf <- unique(demand[c("demand_level", "probability")])
checks <- data.frame(
  check = c("unique support", "nonnegative", "not above one", "total equals one"),
  passes = c(!anyDuplicated(pmf$demand_level), all(pmf$probability >= 0),
             all(pmf$probability <= 1), abs(sum(pmf$probability) - 1) < 1e-10)
)
print(pmf)
print(checks)
stopifnot(all(checks$passes))
print("No probabilities were normalized or changed silently.")

# ----

# PMF and CDF
# Objective: Build PMF and CDF tables and answer exact, cumulative, and tail queries.

# Module 6 R Note 3: PMF and CDF
demand <- read.csv("module6_custom_discrete_demand_synthetic.csv")
pmf <- unique(demand[c("demand_level", "probability")])
pmf <- pmf[order(pmf$demand_level), c("demand_level", "probability")]
pmf$cdf <- cumsum(pmf$probability)
print(pmf)
threshold <- 55
print(c(exact = sum(pmf$probability[pmf$demand_level == threshold]),
        at_most = sum(pmf$probability[pmf$demand_level <= threshold]),
        fewer_than = sum(pmf$probability[pmf$demand_level < threshold]),
        more_than = sum(pmf$probability[pmf$demand_level > threshold])))
barplot(pmf$probability, names.arg = pmf$demand_level, main = "Demand PMF", ylab = "Probability")
plot(pmf$demand_level, pmf$cdf, type = "s", main = "Demand CDF", xlab = "x", ylab = "F(x)")

# ----

# Expected Value and Variance
# Objective: Calculate center, risk, and expected profit from a custom demand PMF.

# Module 6 R Note 4: Expected value and variance
demand <- read.csv("module6_custom_discrete_demand_synthetic.csv")
pmf <- unique(demand[c("demand_level", "probability", "contribution_margin", "shortage_cost", "overage_cost")])
stopifnot(abs(sum(pmf$probability) - 1) < 1e-10)
mu <- sum(pmf$demand_level * pmf$probability)
variance <- sum((pmf$demand_level - mu)^2 * pmf$probability)
order_quantity <- 55
profit <- pmf$demand_level * pmf$contribution_margin -
  pmax(pmf$demand_level - order_quantity, 0) * pmf$shortage_cost -
  pmax(order_quantity - pmf$demand_level, 0) * pmf$overage_cost
print(data.frame(mean = mu, variance = variance, sd = sqrt(variance),
                 expected_profit = sum(profit * pmf$probability)))
print("Expected value is a long-run average, not a guaranteed daily demand.")

# ----

# Discrete Uniform Distribution
# Objective: Compare simulated and theoretical results for equally likely integers.

# Module 6 R Note 5: Discrete uniform
set.seed(62026)
a <- 101
b <- 110
simulated <- sample(a:b, 2000, replace = TRUE)
theoretical_mean <- (a + b) / 2
theoretical_variance <- ((b - a + 1)^2 - 1) / 12
print(data.frame(theoretical_mean, simulated_mean = mean(simulated),
                 theoretical_variance, simulated_variance = var(simulated)))
print(table(simulated))
barplot(table(simulated), main = "Uniform audit draws", xlab = "Invoice number", ylab = "Frequency")

# ----

# Bernoulli and Binomial Models
# Objective: Simulate Bernoulli trials and validate a binomial probability model.

# Module 6 R Note 6: Bernoulli and binomial
services <- read.csv("module6_service_lateness_synthetic.csv")
p <- mean(services$service_late)
n <- 20
x <- 0:n
pmf <- dbinom(x, size = n, prob = p)
print(data.frame(n, p, mean = n * p, sd = sqrt(n * p * (1 - p))))
print(data.frame(x, pmf))
print(sum(pmf))
set.seed(62026)
print(table(rbinom(1000, size = n, prob = p)))
print("The fixed n, constant p, independence, and success-count assumptions require process evidence.")

# ----

# Binomial Probabilities and Shape
# Objective: Translate exact, cumulative, interval, and upper-tail statements.

# Module 6 R Note 7: Binomial probabilities
n <- 24
p <- 0.18
x <- 0:n
results <- data.frame(
  statement = c("exactly 5", "at most 5", "fewer than 5", "more than 5", "at least 5", "between 3 and 7"),
  probability = c(dbinom(5, n, p), pbinom(5, n, p), pbinom(4, n, p),
                  pbinom(5, n, p, lower.tail = FALSE),
                  pbinom(4, n, p, lower.tail = FALSE),
                  pbinom(7, n, p) - pbinom(2, n, p))
)
print(results)
print(data.frame(mean = n * p, variance = n * p * (1 - p), sd = sqrt(n * p * (1 - p))))
plot(x, dbinom(x, n, p), type = "h", lwd = 3, main = "Binomial PMF", xlab = "Late services", ylab = "Probability")
plot(x, pbinom(x, n, p), type = "s", main = "Binomial CDF", xlab = "Late services", ylab = "F(x)")

# ----

# Poisson Counts and Rate Scaling
# Objective: Model arrivals, scale rates to matching units, and compare observed mean and variance.

# Module 6 R Note 8: Poisson
calls <- read.csv("module6_call_arrivals_synthetic.csv")
base_minutes <- 30
lambda <- mean(calls$arrivals[calls$duration_minutes == base_minutes])
target_minutes <- 60
scaled_lambda <- lambda * target_minutes / base_minutes
x <- 0:max(25, ceiling(scaled_lambda + 4 * sqrt(scaled_lambda)))
print(data.frame(observed_mean = mean(calls$arrivals), observed_variance = var(calls$arrivals),
                 base_lambda = lambda, target_minutes, scaled_lambda))
print(c(exactly_10 = dpois(10, scaled_lambda), at_most_10 = ppois(10, scaled_lambda),
        more_than_10 = ppois(10, scaled_lambda, lower.tail = FALSE)))
plot(x, dpois(x, scaled_lambda), type = "h", main = "Poisson PMF", xlab = "Arrivals", ylab = "Probability")
print("Mean equal to variance is a model feature, not automatic proof of a Poisson process.")

# ----

# Poisson Approximation
# Objective: Compare exact binomial and approximate Poisson probabilities.

# Module 6 R Note 9: Poisson approximation to binomial
n <- 120
p <- 0.025
lambda <- n * p
x <- 0:12
exact <- dbinom(x, n, p)
approximate <- dpois(x, lambda)
comparison <- data.frame(x, exact, approximate,
                         absolute_error = abs(exact - approximate),
                         relative_error = ifelse(exact > 0, abs(exact - approximate) / exact, NA))
print(data.frame(n, p, lambda, np = n * p, n_one_minus_p = n * (1 - p)))
print(comparison)
matplot(x, cbind(exact, approximate), type = "b", pch = c(16, 1),
        main = "Exact binomial vs approximate Poisson", xlab = "x", ylab = "Probability")
print("Every Poisson result in this comparison is approximate.")

# ----

# Hypergeometric Sampling
# Objective: Map finite-population parameters and calculate probabilities without replacement.

# Module 6 R Note 10: Hypergeometric
N <- 80
K <- 12
n <- 10
lower <- max(0, n - (N - K))
upper <- min(n, K)
x <- lower:upper
exact <- dhyper(x, m = K, n = N - K, k = n)
print(data.frame(N, K, sample_size = n, support_lower = lower, support_upper = upper))
print(data.frame(x, probability = exact, cdf = phyper(x, m = K, n = N - K, k = n)))
print(sum(exact))
set.seed(62026)
print(table(rhyper(1000, m = K, n = N - K, k = n)))
plot(x, exact, type = "h", lwd = 3, main = "Hypergeometric PMF", xlab = "Defectives sampled", ylab = "Probability")

# ----

# Hypergeometric Approximation
# Objective: Compare exact sampling without replacement with a labelled binomial approximation.

# Module 6 R Note 11: Binomial approximation to hypergeometric
N <- 500
K <- 40
n <- 20
p <- K / N
x <- 0:min(n, K)
exact <- dhyper(x, m = K, n = N - K, k = n)
approximate <- dbinom(x, size = n, prob = p)
comparison <- data.frame(x, exact, approximate,
                         absolute_error = abs(exact - approximate))
print(data.frame(N, K, n, p, sampling_fraction = n / N))
print(comparison)
matplot(x, cbind(exact, approximate), type = "b", pch = c(16, 1),
        main = "Exact hypergeometric vs approximate binomial", xlab = "x", ylab = "Probability")
print("The binomial values are approximate because draws are actually without replacement.")

# ----

# Geometric Waiting Time
# Objective: Use the course trials-until-success convention with base R functions.

# Module 6 R Note 12: Geometric waiting time
p <- 0.22
x <- 1:15
course_pmf <- dgeom(x - 1, prob = p)
course_cdf <- pgeom(x - 1, prob = p)
print(data.frame(course_trials = x, probability = course_pmf, cdf = course_cdf))
print(data.frame(mean_trials = 1 / p, variance_trials = (1 - p) / p^2))
set.seed(62026)
simulated_trials <- rgeom(2000, prob = p) + 1
print(summary(simulated_trials))
plot(x, course_pmf, type = "h", lwd = 3, main = "Trials until first response", xlab = "Course X", ylab = "Probability")
print("Base R counts failures before success; course X counts trials, so use x - 1 and add 1 to simulations.")

# ----

# Transformations, Sums, and Covariance
# Objective: Calculate transformed and aggregate risk under explicit dependence assumptions.

# Module 6 R Note 13: Transformations and covariance
demand <- read.csv("module6_regional_demand_covariance_synthetic.csv")
X <- demand$region_a_demand
Y <- demand$region_b_demand
a <- 12
b <- 450
transformed <- a * X + b
cov_xy <- cov(X, Y)
aggregate_variance <- var(X) + var(Y) + 2 * cov_xy
stopifnot(aggregate_variance >= -1e-10)
print(data.frame(original_mean = mean(X), transformed_mean = mean(transformed),
                 expected_transformed_mean = a * mean(X) + b,
                 transformed_sd = sd(transformed), expected_transformed_sd = abs(a) * sd(X)))
print(data.frame(aggregate_mean = mean(X) + mean(Y),
                 independence_variance = var(X) + var(Y),
                 covariance = cov_xy, aggregate_variance = max(aggregate_variance, 0)))
print("Means always add. Variances require the covariance or a justified independence assumption.")

# ----

# Service Operations Risk Capstone
# Objective: Integrate discrete models into an auditable operating-plan comparison.

# Module 6 R Note 14: Service operations risk capstone
operations <- read.csv("module6_service_risk_capstone.csv")
print(head(operations))
print(str(operations))
print(colSums(is.na(operations)))

arrival_pmf <- prop.table(table(operations$customer_arrivals))
arrival_x <- as.numeric(names(arrival_pmf))
print(data.frame(x = arrival_x, pmf = as.numeric(arrival_pmf),
                 cdf = cumsum(as.numeric(arrival_pmf))))
print(data.frame(expected_arrivals = weighted.mean(arrival_x, arrival_pmf),
                 arrival_sd = sqrt(weighted.mean((arrival_x - weighted.mean(arrival_x, arrival_pmf))^2, arrival_pmf))))

n_services <- 30
p_late <- sum(operations$late_services) / sum(operations$services_completed)
lambda <- mean(operations$customer_arrivals)
print(c(binomial_at_least_6_late = pbinom(5, n_services, p_late, lower.tail = FALSE),
        poisson_more_than_18_arrivals = ppois(18, lambda, lower.tail = FALSE)))

lot <- head(operations, 1)
print(dhyper(2, m = lot$defective_items, n = lot$items_in_lot - lot$defective_items,
             k = lot$audit_sample_size))
conversion_p <- mean(1 / operations$first_conversion_contact)
print(c(geometric_within_5 = pgeom(5 - 1, conversion_p),
        expected_contacts = 1 / conversion_p))

cost <- operations$fixed_staffing_cost +
  operations$variable_service_cost * operations$services_completed +
  operations$shortage_penalty * pmax(operations$customer_arrivals - operations$services_completed, 0)
aggregate_demand <- operations$region_a_demand + operations$region_b_demand
print(data.frame(mean_cost = mean(cost), sd_cost = sd(cost),
                 mean_regional_demand = mean(aggregate_demand),
                 covariance = cov(operations$region_a_demand, operations$region_b_demand)))
par_old <- par(mfrow = c(1, 3))
hist(operations$customer_arrivals, main = "Arrivals", xlab = "Customers")
barplot(arrival_pmf, main = "Empirical arrival PMF", ylab = "Probability")
plot(operations$region_a_demand, operations$region_b_demand,
     main = "Regional demand", xlab = "Region A", ylab = "Region B")
par(par_old)
print("Decision scaffold: compare plans, disclose assumptions and approximations, then write your own recommendation.")

# Original synthetic data generated for STATLAB Academy. No textbook data used.
# Original educational material created for STATLAB Academy. Textbooks may be used only as curriculum references.