# COSC/MATH 201 - Module 8.3: Empirical Models
# Fit functions to NIST data with lm(), check residuals, and compare models.
# Save pontius.csv, danwood.csv, and boxbod.csv in the same folder as this script.

pontius <- read.csv("pontius.csv")
head(pontius)
nrow(pontius)

plot(pontius$load_thousands, pontius$deflection, pch = 16, col = "steelblue",
     xlab = "Load (thousands)", ylab = "Deflection",
     main = "Load cell calibration (NIST Pontius data)")

line_fit <- lm(deflection ~ load_thousands, data = pontius)
coef(line_fit)

plot(pontius$load_thousands, pontius$deflection, pch = 16, col = "steelblue",
     xlab = "Load (thousands)", ylab = "Deflection",
     main = "Least-squares line")
load_grid <- seq(0, 3000, by = 10)
lines(load_grid, predict(line_fit, data.frame(load_thousands = load_grid)),
      col = "firebrick", lwd = 2)

predict(line_fit, data.frame(load_thousands = c(1000, 10000)))

plot(pontius$load_thousands, resid(line_fit), pch = 16, col = "steelblue",
     xlab = "Load (thousands)", ylab = "Residual",
     main = "Residuals from the straight line")
abline(h = 0, lty = 2)

quad_fit <- lm(deflection ~ load_thousands + I(load_thousands^2), data = pontius)
coef(quad_fit)

plot(pontius$load_thousands, resid(quad_fit), pch = 16, col = "darkgreen",
     xlab = "Load (thousands)", ylab = "Residual",
     main = "Residuals from the quadratic")
abline(h = 0, lty = 2)

danwood <- read.csv("danwood.csv")
danwood

plot(danwood$temperature, danwood$energy, pch = 16, cex = 1.3,
     col = "darkorange", xlab = "Temperature (1000 K)",
     ylab = "Energy radiated", main = "Filament data (NIST DanWood)")
lines(range(danwood$temperature), range(danwood$energy), lty = 2, col = "gray50")

powers <- seq(1, 5, by = 0.5)
ssr <- numeric(length(powers))

for (k in seq_along(powers)) {
  z <- danwood$temperature^powers[k]
  fit <- lm(danwood$energy ~ z)
  ssr[k] <- sum(resid(fit)^2)
}

data.frame(power = powers, sum_sq_residuals = round(ssr, 5))

danwood$temp_power <- danwood$temperature^3.5
power_fit <- lm(energy ~ temp_power, data = danwood)
coef(power_fit)

temp_grid <- seq(1.25, 1.75, by = 0.005)
plot(danwood$temperature, danwood$energy, pch = 16, cex = 1.3,
     col = "darkorange", xlab = "Temperature (1000 K)",
     ylab = "Energy radiated", main = "Empirical model f(x) = -0.393 + 0.988 x^3.5")
lines(temp_grid, -0.393 + 0.988 * temp_grid^3.5, col = "firebrick", lwd = 2)

boxbod <- read.csv("boxbod.csv")
boxbod

quad_bod <- lm(bod ~ days + I(days^2), data = boxbod)
poly5_bod <- lm(bod ~ days + I(days^2) + I(days^3) + I(days^4) + I(days^5),
                data = boxbod)
round(resid(poly5_bod), 6)

day_grid <- seq(0.5, 11, by = 0.05)
plot(boxbod$days, boxbod$bod, pch = 16, cex = 1.3, col = "purple",
     xlim = c(0.5, 11), ylim = c(0, 400),
     xlab = "Incubation time (days)", ylab = "BOD (mg/L)",
     main = "Two empirical models for BOD")
lines(day_grid, predict(quad_bod, data.frame(days = day_grid)),
      col = "steelblue", lwd = 2)
lines(day_grid, predict(poly5_bod, data.frame(days = day_grid)),
      col = "firebrick", lwd = 2)
legend("topleft", c("Quadratic", "Degree 5"), col = c("steelblue", "firebrick"),
       lwd = 2, bty = "n")

new_days <- data.frame(days = c(8, 12, 20))
data.frame(days = new_days$days,
           quadratic = round(predict(quad_bod, new_days)),
           degree_5 = round(predict(poly5_bod, new_days)))