
---
title: "Lab: Evaluating Big O Notation for Basic Statistical Methods"
author: "Instructor"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
library(microbenchmark)
library(ggplot2)
```

## 1. Mean Calculation

```{r mean-calculation}
n_values <- c(1e4, 1e5, 1e6, 1e7)

mean_times <- numeric(length(n_values))

for (i in seq_along(n_values)) {
  n <- n_values[i]
  x <- runif(n)
  timing <- microbenchmark(mean_calc = mean(x), times = 5)
  mean_times[i] <- median(timing$time) #/ 1e9
  rm(x)
  gc()
}

mean_df <- data.frame(n = n_values, time = mean_times)

ggplot(mean_df, aes(x = n, y = time)) +
  geom_point() +
  geom_line() +
  scale_x_log10() + scale_y_log10() +
  labs(title = "Mean Calculation Execution Time", x = "Dataset Size (n)", y = "Execution Time (seconds)")

mean_lm <- lm(log10(time) ~ log10(n), data = mean_df)
summary(mean_lm)
```

## 2. Variance Calculation

```{r variance-calculation}
variance_times <- numeric(length(n_values))

for (i in seq_along(n_values)) {
  n <- n_values[i]
  x <- runif(n)
  timing <- microbenchmark(var_calc = var(x), times = 5)
  variance_times[i] <- median(timing$time) / 1e9
  rm(x)
  gc()
}

variance_df <- data.frame(n = n_values, time = variance_times)

ggplot(variance_df, aes(x = n, y = time)) +
  geom_point() +
  geom_line() +
  scale_x_log10() + scale_y_log10() +
  labs(title = "Variance Calculation Execution Time", x = "Dataset Size (n)", y = "Execution Time (seconds)")

variance_lm <- lm(log10(time) ~ log10(n), data = variance_df)
summary(variance_lm)
```

## 3. T-Test Statistic for Difference of Means

```{r ttest-calculation}
ttest_times <- numeric(length(n_values))

for (i in seq_along(n_values)) {
  n_total <- n_values[i]
  n <- n_total / 2
  x1 <- runif(n)
  x2 <- runif(n)
  timing <- microbenchmark(
    t_stat_calc = {
      mean1 <- mean(x1)
      mean2 <- mean(x2)
      var1 <- var(x1)
      var2 <- var(x2)
      t_stat <- (mean1 - mean2) / sqrt(var1 / n + var2 / n)
    },
    times = 5
  )
  ttest_times[i] <- median(timing$time) / 1e9
  rm(x1, x2)
  gc()
}

ttest_df <- data.frame(n = n_values, time = ttest_times)

ggplot(ttest_df, aes(x = n, y = time)) +
  geom_point() +
  geom_line() +
  scale_x_log10() + scale_y_log10() +
  labs(title = "T-Test Statistic Execution Time", x = "Total Sample Size (n)", y = "Execution Time (seconds)")

ttest_lm <- lm(log10(time) ~ log10(n), data = ttest_df)
summary(ttest_lm)
```

## 4. Correlation Coefficient

```{r correlation-calculation}
correlation_times <- numeric(length(n_values))

for (i in seq_along(n_values)) {
  n <- n_values[i]
  x <- runif(n)
  y <- runif(n)
  timing <- microbenchmark(corr_calc = cor(x, y), times = 5)
  correlation_times[i] <- median(timing$time) / 1e9
  rm(x, y)
  gc()
}

correlation_df <- data.frame(n = n_values, time = correlation_times)

ggplot(correlation_df, aes(x = n, y = time)) +
  geom_point() +
  geom_line() +
  scale_x_log10() + scale_y_log10() +
  labs(title = "Correlation Coefficient Execution Time", x = "Dataset Size (n)", y = "Execution Time (seconds)")

correlation_lm <- lm(log10(time) ~ log10(n), data = correlation_df)
summary(correlation_lm)
```

## 5. Linear Regression (for n)

```{r regression-calculation}
p <- 5
regression_times <- numeric(length(n_values))

for (i in seq_along(n_values)) {
  n <- n_values[i]
  X <- matrix(runif(n * p), nrow = n, ncol = p)
  y <- runif(n)
  timing <- microbenchmark(
    lm_calc = {
      XtX <- crossprod(X)
      XtX_inv <- solve(XtX)
      Xty <- crossprod(X, y)
      beta_hat <- XtX_inv %*% Xty
    },
    times = 5
  )
  regression_times[i] <- median(timing$time) / 1e9
  rm(X, y, XtX, XtX_inv, Xty, beta_hat)
  gc()
}

regression_df <- data.frame(n = n_values, time = regression_times)

ggplot(regression_df, aes(x = n, y = time)) +
  geom_point() +
  geom_line() +
  scale_x_log10() + scale_y_log10() +
  labs(title = "Linear Regression Execution Time", x = "Dataset Size (n)", y = "Execution Time (seconds)")

regression_lm <- lm(log10(time) ~ log10(n), data = regression_df)
summary(regression_lm)
```
## 5. Linear Regression (for p)

```{r}
library(microbenchmark)
library(ggplot2)

# Define predictor count values
p_vals <- c(5, 10, 15, 20,25,30,35,40)
regression_times <- numeric(length(p_vals))  # Use p_vals to initialize length

for (i in seq_along(p_vals)) {
  n <- 1e4  # Fixed number of samples
  p <- p_vals[i]  # Set current predictor count
  X <- matrix(runif(n * p), nrow = n, ncol = p)
  y <- runif(n)
  
  timing <- microbenchmark(
    lm_calc = {
      XtX <- crossprod(X)
      XtX_inv <- solve(XtX)
      Xty <- crossprod(X, y)
      beta_hat <- XtX_inv %*% Xty
    },
    times = 5
  )
  
  regression_times[i] <- median(timing$time) / 1e9  # Convert nanoseconds to seconds
  rm(X, y, XtX, XtX_inv, Xty, beta_hat)  # Clear memory
  gc()
}

# Create data frame
regression_df <- data.frame(p = p_vals, time = regression_times)

# Plot
ggplot(regression_df, aes(x = p, y = time)) +
  geom_point() +
  geom_line() +
  labs(title = "Linear Regression Execution Time", x = "Predictor Count (p)", y = "Execution Time (seconds)")

# Fit a linear model to check relationship (use `p` instead of `p_vals` in the data)
regression_lm <- lm(time ~ p, data = regression_df)
summary(regression_lm)

```

## 6. Principal Component Analysis (PCA)

```{r pca-calculation}
n_values_pca <- c(100, 200, 300,400,500,600)
pca_times <- numeric(length(n_values_pca))

for (i in seq_along(n_values_pca)) {
  n <- n_values_pca[i]
  p <- n
  X <- matrix(rnorm(n * p), nrow = n, ncol = p)
  timing <- microbenchmark(
    pca_calc = {
      X_centered <- scale(X, center = TRUE, scale = FALSE)
      C <- cov(X_centered)
      eig <- eigen(C, symmetric = TRUE)
      k <- 5
      principal_components <- X_centered %*% eig$vectors[, 1:k]
    },
    times = 3
  )
  pca_times[i] <- median(timing$time) / 1e9
  rm(X, X_centered, C, eig, principal_components)
  gc()
}

pca_df <- data.frame(n = n_values_pca, time = pca_times)

ggplot(pca_df, aes(x = n, y = time)) +
  geom_point() +
  geom_line() +
  labs(title = "PCA Execution Time", x = "Dataset Size (n = p)", y = "Execution Time (seconds)")

pca_lm <- lm(log10(time) ~ log10(n), data = pca_df)
summary(pca_lm)
```
