# Continuous Variables, PDFs, and CDFs
# Objective: Connect density, cumulative area, and zero point probability.

# Module 7 R Note 1: Continuous foundations
service <- read.csv("module7_service_times_synthetic.csv")
x <- service$service_minutes
print(summary(x))
print(c(sample_mean = mean(x), sample_sd = sd(x)))
grid <- seq(min(x), max(x), length.out = 200)
density_estimate <- density(x, from = min(grid), to = max(grid))
plot(density_estimate, main = "Synthetic service-time density", xlab = "Minutes")
abline(v = c(20, 30), col = "#d95d39", lty = 2)
print("For a continuous model P(X = 25) = 0; interval probability is area.")

# ----

# Continuous Uniform Distribution
# Objective: Calculate bounded interval probabilities and compare theory with simulation.

# Module 7 R Note 2: Continuous uniform
a <- 18
b <- 42
lower <- 24
upper <- 33
probability <- max(0, min(upper, b) - max(lower, a)) / (b - a)
print(data.frame(a, b, mean = (a + b) / 2, sd = (b - a) / sqrt(12), probability))
set.seed(72026)
draws <- runif(5000, a, b)
print(c(simulated = mean(draws >= lower & draws <= upper), theoretical = probability))
hist(draws, probability = TRUE, main = "Uniform completion time", xlab = "Days")
abline(h = 1 / (b - a), col = "#087f8c", lwd = 2)

# ----

# Normal Shape and Parameters
# Objective: Explore how mean and standard deviation control a normal curve.

# Module 7 R Note 3: Normal shape
packages <- read.csv("module7_package_weights_synthetic.csv")
x <- packages$weight_grams
mu <- mean(x)
sigma <- sd(x)
grid <- seq(mu - 4 * sigma, mu + 4 * sigma, length.out = 300)
print(data.frame(mean = mu, sd = sigma, within_1_sd = mean(abs(x - mu) <= sigma), within_2_sd = mean(abs(x - mu) <= 2 * sigma)))
hist(x, probability = TRUE, breaks = 14, main = "Package weights", xlab = "Grams")
lines(grid, dnorm(grid, mu, sigma), col = "#087f8c", lwd = 3)
print("Symmetry and bell shape must be checked; they are not guaranteed by a mean and SD.")

# ----

# Normal Probability and z Scores
# Objective: Calculate left, right, and interval probabilities without tail errors.

# Module 7 R Note 4: Normal probability
mu <- 500
sigma <- 12
lower <- 485
upper <- 520
results <- c(left = pnorm(lower, mu, sigma), right = pnorm(upper, mu, sigma, lower.tail = FALSE), between = pnorm(upper, mu, sigma) - pnorm(lower, mu, sigma))
print(results)
print(c(z_lower = (lower - mu) / sigma, z_upper = (upper - mu) / sigma))
curve(dnorm(x, mu, sigma), from = mu - 4 * sigma, to = mu + 4 * sigma, main = "Normal package model", xlab = "Grams")
print("For continuous X, strict and inclusive interval endpoints have the same probability.")

# ----

# Inverse Normal and Percentiles
# Objective: Find thresholds from cumulative probabilities and central coverage.

# Module 7 R Note 5: Inverse normal
mu <- 500
sigma <- 12
left_area <- 0.95
upper_threshold <- qnorm(left_area, mu, sigma)
central <- 0.90
tail <- (1 - central) / 2
central_limits <- qnorm(c(tail, 1 - tail), mu, sigma)
print(data.frame(left_area, upper_threshold))
print(data.frame(central_coverage = central, lower = central_limits[1], upper = central_limits[2]))
print(pnorm(upper_threshold, mu, sigma))
print("qnorm expects a left-tail cumulative probability.")

# ----

# Normal Approximations
# Objective: Compare exact count probabilities with continuity-corrected normal approximations.

# Module 7 R Note 6: Normal approximations
n <- 160
p <- 0.12
cutoff <- 25
mu_bin <- n * p
sd_bin <- sqrt(n * p * (1 - p))
exact_binomial <- pbinom(cutoff - 1, n, p, lower.tail = FALSE)
approx_binomial <- pnorm(cutoff - 0.5, mu_bin, sd_bin, lower.tail = FALSE)
lambda <- 36
poisson_cutoff <- 44
exact_poisson <- ppois(poisson_cutoff, lambda)
approx_poisson <- pnorm(poisson_cutoff + 0.5, lambda, sqrt(lambda))
print(data.frame(model = c("binomial", "Poisson"), exact = c(exact_binomial, exact_poisson), approximate = c(approx_binomial, approx_poisson), absolute_error = abs(c(exact_binomial, exact_poisson) - c(approx_binomial, approx_poisson))))
print(c(np = n * p, n_one_minus_p = n * (1 - p), lambda = lambda))
print("Every normal approximation is labelled approximate and uses a 0.5 continuity correction.")

# ----

# Exponential Waiting Times
# Objective: Calculate waiting-time probabilities under a constant-rate process.

# Module 7 R Note 7: Exponential waiting time
arrivals <- read.csv("module7_customer_arrivals_synthetic.csv")
lambda <- 1 / mean(arrivals$minutes_until_next_arrival)
lower <- 2
upper <- 7
results <- c(within_upper = pexp(upper, lambda), longer_than_upper = pexp(upper, lambda, lower.tail = FALSE), between = pexp(upper, lambda) - pexp(lower, lambda))
print(data.frame(rate_per_minute = lambda, mean_wait = 1 / lambda))
print(results)
curve(dexp(x, lambda), from = 0, to = qexp(0.995, lambda), main = "Waiting-time density", xlab = "Minutes")
print("The constant event-rate and independence assumptions require process evidence.")

# ----

# Inverse Exponential, MTBF, and Warranty
# Objective: Find service levels and warranty thresholds from a target probability.

# Module 7 R Note 8: Inverse exponential and warranty
units <- read.csv("module7_warranty_lifetimes_synthetic.csv")
mtbf <- mean(units$lifetime_hours)
rate <- 1 / mtbf
failure_target <- 0.10
warranty <- qexp(failure_target, rate)
service_level <- 0.95
response_threshold <- qexp(service_level, rate)
print(data.frame(mtbf, rate, warranty_at_10_percent_failure = warranty, threshold_at_95_percent = response_threshold))
print(c(check_warranty = pexp(warranty, rate), survival_at_warranty = pexp(warranty, rate, lower.tail = FALSE)))
print("A constant hazard is a strong assumption; inspect aging, maintenance, and censoring.")

# ----

# Triangular What-if Model
# Objective: Use bounded minimum, most-likely, and maximum values for transparent scenarios.

# Module 7 R Note 9: Triangular what-if model
rtri <- function(n, a, b, c) {
  u <- runif(n)
  split <- (b - a) / (c - a)
  ifelse(u < split, a + sqrt(u * (b - a) * (c - a)), c - sqrt((1 - u) * (c - b) * (c - a)))
}
set.seed(72026)
a <- 20; b <- 45; c <- 90
draws <- rtri(10000, a, b, c)
print(data.frame(theoretical_mean = (a + b + c) / 3, simulated_mean = mean(draws), simulated_sd = sd(draws), probability_under_60 = mean(draws <= 60)))
hist(draws, probability = TRUE, breaks = 30, main = "Triangular project duration", xlab = "Days")
print("The parameters are stakeholder scenarios, not fitted facts.")

# ----

# Continuous Model Selection
# Objective: Compare candidate models using support, shape, process, and decision evidence.

# Module 7 R Note 10: Model selection
service <- read.csv("module7_service_times_synthetic.csv")
x <- service$service_minutes
diagnostics <- data.frame(n = length(x), minimum = min(x), mean = mean(x), median = median(x), maximum = max(x), sd = sd(x), cv = sd(x) / mean(x))
print(diagnostics)
print(c(normal_tail_above_40 = pnorm(40, mean(x), sd(x), lower.tail = FALSE), empirical_tail_above_40 = mean(x > 40)))
par_old <- par(mfrow = c(1, 2))
hist(x, probability = TRUE, main = "Observed shape", xlab = "Minutes")
qqnorm(x); qqline(x, col = "#d95d39")
par(par_old)
print("Select a model from process logic and diagnostics; do not choose from appearance alone.")

# ----

# Northstar Operations Risk Lab
# Objective: Integrate continuous models into an auditable service, quality, and warranty recommendation.

# Module 7 R Note 11: Northstar Operations Risk Lab
northstar <- read.csv("module7_northstar_capstone.csv")
print(head(northstar))
print(colSums(is.na(northstar)))
service_mu <- mean(northstar$service_minutes)
service_sd <- sd(northstar$service_minutes)
wait_rate <- 1 / mean(northstar$minutes_until_next_arrival)
life_rate <- 1 / mean(northstar$lifetime_hours)
quality_n <- 120
quality_p <- mean(northstar$defect_flag)
exact_quality <- pbinom(18, quality_n, quality_p, lower.tail = FALSE)
approx_quality <- pnorm(18.5, quality_n * quality_p, sqrt(quality_n * quality_p * (1 - quality_p)), lower.tail = FALSE)
project_draws <- northstar$project_duration_days
summary_table <- data.frame(metric = c("Service above 35 minutes", "Wait above 8 minutes", "Warranty failure within 1000 hours", "Exact quality tail", "Approximate quality tail", "Project P90"), value = c(pnorm(35, service_mu, service_sd, lower.tail = FALSE), pexp(8, wait_rate, lower.tail = FALSE), pexp(1000, life_rate), exact_quality, approx_quality, unname(quantile(project_draws, 0.90))))
print(summary_table)
print(data.frame(service_mu, service_sd, wait_rate, life_rate, quality_p, approximation_error = abs(exact_quality - approx_quality)))
par_old <- par(mfrow = c(1, 3))
hist(northstar$service_minutes, main = "Service", xlab = "Minutes")
hist(northstar$minutes_until_next_arrival, main = "Waiting", xlab = "Minutes")
hist(project_draws, main = "Project scenarios", xlab = "Days")
par(par_old)
print("Student decision: compare thresholds, name assumptions, quantify approximation error, and write an original recommendation.")

# Original synthetic data generated for STATLAB Academy. No textbook data used.
# Educational material for STATLAB Academy. Verify assumptions, units, tail direction, and rounding before making a decision.