set.seed(1)
library(readxl)
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(gridExtra)

## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine

source('helper.R')
# read in the data
df <- read_excel('kidney_data.xlsx')

dage <- df$dage


# extract population parameters (N, mu, sigma_sq)
N <- nrow(df)

mu <- mean(dage)

sigma_sq <- sum((dage - mu)^2) / N



# ---------------------------------------
# problem 2


n = 5

n_sims <- 100

simulation_dataset <- matrix(nrow = 5, ncol = 100)

for(i in 1:n_sims){
  sample_i <- sample(dage, n, replace = TRUE)
  simulation_dataset[,i] <- sample_i
}

# display only the first 4 simulations

simulation_display <- simulation_dataset[,1:4]
simulation_display

##      [,1] [,2] [,3] [,4]
## [1,]   48   20   26   39
## [2,]   29   40   25   24
## [3,]   32   24   32   39
## [4,]   48   38   69   24
## [5,]   42   42   37   19

# problem 2b

x_bars <- apply(simulation_dataset, 2, mean) 

expectation_x_bar <- mean(x_bars)
expectation_x_bar

## [1] 33.822

pop_mean <- mean(df$dage)
pop_mean

## [1] 33.63194

# problem 2c

empirical_sampling_variance <- sd(x_bars)
empirical_sampling_variance

## [1] 5.104476

# problem 2d

theoretical_sampling_variance <- sqrt(sigma_sq) / sqrt(n)
theoretical_sampling_variance

## [1] 5.108475

# problem 2e

sample_variances <- apply(simulation_dataset, 2, var)

expectation_sample_variances <- mean(sample_variances)
expectation_sample_variances

## [1] 130.259

sigma_sq

## [1] 130.4826

#Problem 2g
empirical_s_sq_std_err <- sd(sample_variances) 
empirical_s_sq_std_err

## [1] 84.69829

theoretical_s_sq_std_err <- sqrt(2 * sigma_sq^2 / (n-1)) 
theoretical_s_sq_std_err

## [1] 92.26512

# -----------------------------------------
# Problem 3

n_5_results  <- generate_simulations(dage, n = 5,  n_sims = 100)
n_5_results

## x_bar expectation     x_bar std err   s^2 expectation       s^2 std err 
##          33.68400           5.17056         137.44100          92.60296

n_10_results <- generate_simulations(dage, n = 10, n_sims = 100)
n_10_results

## x_bar expectation     x_bar std err   s^2 expectation       s^2 std err 
##         33.971000          3.330614        142.902333         63.761615

n_50_results <- generate_simulations(dage, n = 50, n_sims = 100)
n_50_results

## x_bar expectation     x_bar std err   s^2 expectation       s^2 std err 
##         33.596800          1.398226        131.977494         25.624357

# problem 4

generate_histograms(dage, 100, mu, sigma_sq)

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).

## TableGrob (3 x 2) "arrange": 6 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## 5 5 (3-3,1-1) arrange gtable[layout]
## 6 6 (3-3,2-2) arrange gtable[layout]