knitr::opts_chunk$set(echo = TRUE)
#I have read the honor code and abided by it.
# I strive to uphold the University values of respect, responsibility, discovery, and excellence. On my honor, I pledge that I have neither given nor received unauthorized assistance on this work.
#Question 1A
x <- c(12, NA, 5, 20, 0, 9)
x[x != NA & x >= 9]
## [1] NA NA NA NA
#!= NA returns all NAs because R can't compare anything directly to NA using !=
retained <- x[!is.na(x) & x >= 9]
retained
## [1] 12 20 9
mean(retained)
## [1] 13.66667
#Question 1B
weight <- factor(c("10", "5", "20", "10"),
levels = c("5", "10", "20"))
as.numeric(weight)
## [1] 2 1 3 2
#as.numeric() pulls levels, not the actual values, so it needs to be made a character first
weight_k<- as.numeric(as.character(weight))
weight_k
## [1] 10 5 20 10
weight_g <- weight_k * 1000
weight_g
## [1] 10000 5000 20000 10000
#Question 1C
A <- matrix(1:12, nrow = 4, ncol = 3)
A
## [,1] [,2] [,3]
## [1,] 1 5 9
## [2,] 2 6 10
## [3,] 3 7 11
## [4,] 4 8 12
AxA<- A * A
AxA
## [,1] [,2] [,3]
## [1,] 1 25 81
## [2,] 4 36 100
## [3,] 9 49 121
## [4,] 16 64 144
dim(AxA)
## [1] 4 3
AtA <- t(A) %*% A
AtA
## [,1] [,2] [,3]
## [1,] 30 70 110
## [2,] 70 174 278
## [3,] 110 278 446
dim(AtA)
## [1] 3 3
#In AxA each entry in the matrix is multiplied by itself,which is why it stays a 4x3, but AtA multiplies a 3x4 and 4x3 matrix making a 3x3 matrix
col_means <- apply(A, 2, mean)
col_means
## [1] 2.5 6.5 10.5
#No this does not have an inverse because column 3 depends on columns 1 and 2
#aka the columns are not indepedent
#Question 1D
matrix_info <- list(
data = A,
means = col_means,
crossproduct = AtA)
str(matrix_info["means"])
## List of 1
## $ means: num [1:3] 2.5 6.5 10.5
str(matrix_info[["means"]])
## num [1:3] 2.5 6.5 10.5
#The single brackets return a list, while double brackets pull out the actual object
matrix_info[["crossproduct"]][2, 3]
## [1] 278
#Question 2A
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.3 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)
depots <- read.csv("depots.csv",
colClasses = c(depot_id = "character"),
na.strings = "M")
deliveries <- read.csv("deliveries.csv",
colClasses = c(depot_id = "character"),
na.strings = "M")
deliveries$date <- as.Date(deliveries$date)
dim(depots)
## [1] 5 2
dim(deliveries)
## [1] 180 7
str(deliveries)
## 'data.frame': 180 obs. of 7 variables:
## $ delivery_id: chr "D001" "D002" "D003" "D004" ...
## $ depot_id : chr "04" "04" "03" "01" ...
## $ date : Date, format: "2026-07-20" "2026-09-11" ...
## $ distance_km: num 19 4 24.9 18.3 17.5 12.8 NA 20.6 16.8 19.1 ...
## $ weight_kg : num 17.7 11.6 7.1 3.6 7.2 1.1 11.3 10.4 7.8 2.4 ...
## $ minutes : num 86.6 31.7 86.1 67 49.1 49.2 33 77.9 69.5 70.6 ...
## $ rating : int 4 5 4 2 NA 3 5 1 4 4 ...
#Question 2B
delivery<-deliveries |>
left_join(depots, by = "depot_id")
nrow(delivery) == nrow(deliveries)
## [1] TRUE
!any(is.na(delivery$depot))
## [1] TRUE
unused_depot <- depots |>
filter(!depot_id %in% delivery$depot_id) |>
pull(depot)
unused_depot
## [1] "Central"
delivery$depot <- factor(delivery$depot, levels = depots$depot)
#Question 2C
#I used AI on this one because it was picking up minutes from the lubricate package, so I had figure that out
depot_summary <- delivery |>
group_by(depot) |>
summarise(
n_deliveries = n(),
n_observed_times = sum(!is.na(minutes)),
mean_time = mean(minutes, na.rm = TRUE),
sd_time = sd(minutes, na.rm = TRUE),
.groups = "drop"
) |>
arrange(desc(mean_time))
print(depot_summary)
## # A tibble: 4 × 5
## depot n_deliveries n_observed_times mean_time sd_time
## <fct> <int> <int> <dbl> <dbl>
## 1 South 45 45 79.1 31.8
## 2 West 45 44 74.7 27.9
## 3 North 45 44 74.1 29.6
## 4 East 45 43 66.5 28.8
write.csv(depot_summary, file = "depot_summary.csv", row.names = FALSE)
#The number of deliveries and observe ones differ due to missing data. The observed times only count it when the variable is present.
#Question 2D
delivery_long <- delivery |>
pivot_longer(
cols = c(distance_km, weight_kg, minutes),
names_to = "measurement",
values_to = "value")|>
select(delivery_id, measurement, value)
dim(delivery_long)
## [1] 540 3
delivery_wide <- delivery_long |>
pivot_wider(
names_from = measurement,
values_from = value)
delivery_long_ordered<- delivery |>
select(delivery_id, distance_km, weight_kg, minutes) |>
arrange(delivery_id)
all.equal(delivery_long_ordered, delivery_wide)
## [1] "Attributes: < Component \"class\": Lengths (1, 3) differ (string compare on first 1) >"
## [2] "Attributes: < Component \"class\": 1 string mismatch >"
#Question 2 E
plot(x = delivery$minutes, y = delivery$distance_km)
library(ggplot2)
pl <-ggplot(delivery, aes(x = minutes, y = distance_km)) +
geom_point() +
facet_wrap(~depot)+
geom_smooth(method = "lm", se = FALSE)
pl
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 9 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).
pl_plot<- file.path(tempdir(), "delivery_plot.png")
ggsave(filename = "delivery_plot.png", plot = pl, width = 7, height = 5, units = "in", dpi = 150)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 9 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).
file.exists(pl_plot)
## [1] FALSE
#Yes this relationship looks similar across the depots.
#Question 3A
delivery2<-delivery[!is.na(delivery$minutes) &
!is.na(delivery$distance_km) &
!is.na(delivery$weight_kg), ]
nrow(delivery2)
## [1] 165
set.seed(202603)
train_amo<- sample(1:165, size = 120, replace = FALSE)
train <- delivery2[train_amo, ]
test <- delivery2[-train_amo, ]
nrow(train)
## [1] 120
nrow(test)
## [1] 45
m1<-lm(minutes ~ distance_km + weight_kg, data = train)
summary(m1)
##
## Call:
## lm(formula = minutes ~ distance_km + weight_kg, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.159 -4.054 -0.574 4.615 19.500
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.85638 1.76905 6.137 1.18e-08 ***
## distance_km 3.02914 0.06858 44.173 < 2e-16 ***
## weight_kg 1.09019 0.12843 8.489 7.56e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.002 on 117 degrees of freedom
## Multiple R-squared: 0.9451, Adjusted R-squared: 0.9441
## F-statistic: 1007 on 2 and 117 DF, p-value: < 2.2e-16
#This is below .05 so it is SS. No, this does not establish causation as there are other factors which may impact the results.
#Question 3B
X <- cbind(1, train$distance_km, train$weight_kg)
y <- train$minutes
math<-solve(t(X) %*% X) %*% t(X) %*% y
dim(X)
## [1] 120 3
dim(t(X) %*% X)
## [1] 3 3
all.equal(as.vector(math), as.vector(coef(m1)))
## [1] TRUE
#you cannot use == as the variables are not comparable between each other
#Question 3C
new<- data.frame(distance_km = 15, weight_kg = 5)
predict(m1, newdata = new, interval = "confidence", level = 0.95)
## fit lwr upr
## 1 61.74439 60.09383 63.39495
predict(m1, newdata = new, interval = "prediction", level = 0.95)
## fit lwr upr
## 1 61.74439 47.77933 75.70945
#The prediction interval is wider as it accounts for uncertanty
par(mfrow = c(1, 2))
plot(m1)
#These plots help to show that regression diagnostics have been pass!
#Question 3D
#I used AI to help calculate the RSME. I got confused with the equation, so I asked it to help me with that and making the table with everything.
library(randomForest)
## randomForest 4.7-1.2
## Type rfNews() to see new features/changes/bug fixes.
##
## Attaching package: 'randomForest'
##
## The following object is masked from 'package:dplyr':
##
## combine
##
## The following object is masked from 'package:ggplot2':
##
## margin
set.seed(2026031)
rf<- randomForest(minutes ~ distance_km + weight_kg,
data = train,
ntree = 300,
mtry = 2)
rmse <- function(actual, predicted) {
sqrt(mean((actual - predicted)^2))}
mean_y <- mean(train$minutes)
train_rmse_base <- rmse(train$minutes, mean_y)
test_rmse_base <- rmse(test$minutes, mean_y)
train_rmse_lm <- rmse(train$minutes, predict(m1, newdata = train))
test_rmse_lm <- rmse(test$minutes, predict(m1, newdata = test))
train_rmse_rf <- rmse(train$minutes, predict(rf, newdata = train))
test_rmse_rf <- rmse(test$minutes, predict(rf, newdata = test))
data.frame(
Method = c("Baseline", "Linear Model", "Random Forest"),
`Training RMSE` = c(train_rmse_base, train_rmse_lm, train_rmse_rf),
`Test RMSE` = c(test_rmse_base, test_rmse_lm, test_rmse_rf))
## Method Training.RMSE Test.RMSE
## 1 Baseline 29.503495 29.860557
## 2 Linear Model 6.913960 7.112118
## 3 Random Forest 3.491358 8.821323
#Question 3E
#AI helped me with the bootstrap. I understand the logic of this, but still have not gotten the hang of the actual code so i worked it though with it.
set.seed(2026032)
B <- 2000
n_train <- nrow(train)
boot_coefs <- numeric(B)
for (b in 1:B) {
boot_indices <- sample(1:n_train, size = n_train, replace = TRUE)
boot_sample <- train[boot_indices, ]
boot_fit <- lm(minutes ~ distance_km + weight_kg, data = boot_sample)
boot_coefs[b] <- coef(boot_fit)["distance_km"]}
hist(boot_coefs)
quantile<-quantile(boot_coefs, probs = c(0.025, 0.975))
quantile
## 2.5% 97.5%
## 2.884031 3.173399
confint<-confint(m1, parm = "distance_km", level = 0.95)
confint
## 2.5 % 97.5 %
## distance_km 2.893328 3.164947
#resampling would be incorrect because it breaks down the relationship between the variables
#Question 4A
weight <- factor(c("10", "5", "20", "10"),
levels = c("5", "10", "20"))
as.numeric(weight)
## [1] 2 1 3 2
#as.numeric() pulls levels, not the actual values, so it needs to be made a character first
d_target <- function(x) {
ifelse(x >= 0 & x <= 1, 6 * x * (1 - x), 0)
}
test<- c(-0.1, 0, 0.25, 0.5, 1, 1.1)
d_target(test)
## [1] 0.000 0.000 1.125 1.500 0.000 0.000
#Question 4B
# I used AI to help me write the loop. I understand how to read the code, but have struggled to actually write the code, so it helped me wit that step by step.
r_target <- function(n) {
if (n < 1 || n %% 1 != 0) stop("n >= 1")
draws <- numeric(n)
accepted <- 0
attempts <- 0
while (accepted < n) {
attempts <- attempts + 1
u <- runif(1); v <- runif(1)
if (v <= 4 * u * (1 - u)) {
accepted <- accepted + 1
draws[accepted] <- u}}
list(draws = draws, acceptance = n, acceptance = n / attempts)
}
#Question 4C
set.seed(202604)
res <- r_target(5000)
mean(res$draws)
## [1] 0.5014082
var(res$draws)
## [1] 0.04984302
res$acceptance
## [1] 5000
res$attempts
## NULL
hist(res$draws)
#Question 4D
#AI helped me calculate the MCSE and transpose the equation into code. Like, I understand the equation, but need to be walkthrough how to actually make it into code
p_hat <- mean(res$draws > 0.8)
mcse <- sqrt(p_hat * (1 - p_hat) / 5000)
p<- pbeta(0.8, 2, 2)
q<- qbeta(0.90, 2, 2)
p_hat
## [1] 0.1048
mcse
## [1] 0.004331673
p
## [1] 0.896
q
## [1] 0.8041999
#Question 5A
running_mean_loop <- function(x) {
out <- numeric(length(x))
for (i in seq_along(x)) {
out[i] <- mean(x[seq_len(i)])
}
out
}
running_mean_vectorized <- function(x) {
cumsum(x) / seq_along(x)
}
test2<- c(2, 4, 9, 1)
running_mean_loop(test2)
## [1] 2 3 5 4
running_mean_vectorized(test2)
## [1] 2 3 5 4
running_mean_loop(numeric(0))
## numeric(0)
running_mean_vectorized(numeric(0))
## numeric(0)
set.seed(202605)
x <- rnorm(10000)
loop <- Sys.time()
for (i in 1:5) running_mean_loop(x)
difftime(Sys.time(), loop)
## Time difference of 1.11096 secs
vec <- Sys.time()
for (i in 1:5) running_mean_vectorized(x)
difftime(Sys.time(), vec)
## Time difference of 0.0008971691 secs
#The loop takes ~0.20 seconds and the vectorized version takes ~0.0003 seconds. The vectorization is not always faster as it can cause increased memory making it go slower
#Question 5B
#AI helped me write this code as I was unfamiliar with this package and needed to be walked through it
library(future.apply)
## Loading required package: future
plan(sequential)
means_seq <- unlist(future_lapply(
1:8,
function(i) mean(r_target(1000)),
future.seed = 2026051
))
## Warning in mean.default(r_target(1000)): argument is not numeric or logical:
## returning NA
## Warning in mean.default(r_target(1000)): argument is not numeric or logical:
## returning NA
## Warning in mean.default(r_target(1000)): argument is not numeric or logical:
## returning NA
## Warning in mean.default(r_target(1000)): argument is not numeric or logical:
## returning NA
## Warning in mean.default(r_target(1000)): argument is not numeric or logical:
## returning NA
## Warning in mean.default(r_target(1000)): argument is not numeric or logical:
## returning NA
## Warning in mean.default(r_target(1000)): argument is not numeric or logical:
## returning NA
## Warning in mean.default(r_target(1000)): argument is not numeric or logical:
## returning NA
workers <- min(2, availableCores())
plan(multisession, workers = workers)
means_multi <- unlist(future_lapply(
1:8,
function(i) mean(r_target(1000)),
future.seed = 2026051
))
## Warning in mean.default(r_target(1000)): argument is not numeric or logical:
## returning NA
## Warning in mean.default(r_target(1000)): argument is not numeric or logical:
## returning NA
## Warning in mean.default(r_target(1000)): argument is not numeric or logical:
## returning NA
## Warning in mean.default(r_target(1000)): argument is not numeric or logical:
## returning NA
## Warning in mean.default(r_target(1000)): argument is not numeric or logical:
## returning NA
## Warning in mean.default(r_target(1000)): argument is not numeric or logical:
## returning NA
## Warning in mean.default(r_target(1000)): argument is not numeric or logical:
## returning NA
## Warning in mean.default(r_target(1000)): argument is not numeric or logical:
## returning NA
means_seq
## [1] NA NA NA NA NA NA NA NA
means_multi
## [1] NA NA NA NA NA NA NA NA
sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.7.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/Chicago
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] future.apply_1.21.0 future_1.80.0 randomForest_4.7-1.2
## [4] lubridate_1.9.5 forcats_1.0.1 stringr_1.6.0
## [7] dplyr_1.2.1 purrr_1.2.2 readr_2.2.0
## [10] tidyr_1.3.2 tibble_3.3.1 ggplot2_4.0.3
## [13] tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.10 utf8_1.2.6 generics_0.1.4 stringi_1.8.9
## [5] lattice_0.23-1 listenv_1.1.0 hms_1.1.4 digest_0.6.39
## [9] magrittr_2.0.5 evaluate_1.0.5 grid_4.6.1 timechange_0.4.0
## [13] RColorBrewer_1.1-3 fastmap_1.2.0 jsonlite_2.0.0 Matrix_1.7-6
## [17] mgcv_1.9-4 scales_1.4.0 codetools_0.2-20 textshaping_1.0.5
## [21] jquerylib_0.1.4 cli_3.6.6 rlang_1.3.0 parallelly_1.48.0
## [25] splines_4.6.1 withr_3.0.3 cachem_1.1.0 yaml_2.3.12
## [29] otel_0.2.0 parallel_4.6.1 tools_4.6.1 tzdb_0.5.0
## [33] globals_0.20.0 vctrs_0.7.3 R6_2.6.1 lifecycle_1.0.5
## [37] ragg_1.5.2 pkgconfig_2.0.3 pillar_1.11.1 bslib_0.12.0
## [41] gtable_0.3.6 glue_1.8.1 systemfonts_1.3.2 xfun_0.62
## [45] tidyselect_1.2.1 rstudioapi_0.19.0 knitr_1.52 farver_2.1.2
## [49] htmltools_0.5.9 nlme_3.1-171 rmarkdown_2.32 labeling_0.4.3
## [53] compiler_4.6.1 S7_0.2.2