## ----setup, include=FALSE-----------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)


## ----clt----------------------------------------------------------------------------------------------
##draw n=5 samples from uniform distribution and repeat this 
##experiment 200 times 
sim1=replicate(200,mean(runif(5)))
hist(sim1,xlim=c(0,1),freq = F, main="N=5")
lines(density(sim1),col="red")

##draw n=100 samples from uniform distribution and repeat this 
##experiment 200 times 
sim2=replicate(200,mean(runif(100)))
hist(sim2,freq = F, main="N=100", xlim=c(0, 1))
lines(density(sim2),col="red")



## ----q2-----------------------------------------------------------------------------------------------
library(dslabs)
data(murders)
boxplot(population~region, data=murders)


## -----------------------------------------------------------------------------------------------------
data("ChickWeight")
str(ChickWeight)
ChickWeight[1:5,]


## -----------------------------------------------------------------------------------------------------
library(NHANES)
data(NHANES)
subset_data <- subset(NHANES, Gender == "female" & AgeDecade == " 20-29")
mean_bp <- mean(subset_data$BPSysAve, na.rm = TRUE)  # Remove missing values
sd_bp <- sd(subset_data$BPSysAve, na.rm = TRUE)      # Remove missing values
ref <- c(mean = mean_bp, sd = sd_bp)

print(ref)


## -----------------------------------------------------------------------------------------------------
library(tidyverse)
url <- "https://people.stat.sc.edu/hoyen/STAT540/Data/baby.txt"
babiesraw <- read.table(file = url, header = TRUE)
babies<-filter(babiesraw, smoke!=9)
str(babies)
stripchart(bwt ~ smoke, data = babies, method = "jitter", jitter = 0.2, vertical = TRUE, 
           pch = 1, xlab = "Smoking Status (0 = non-smoker, 1 = smoker)", 
           ylab = "Birth Weight (bwt)", main = "Scatter plots for babies' birth weight")


x1 <- 0.9  
x2 <- 2.1  
y1 <- mean(babies$bwt[babies$smoke == 0], na.rm = TRUE) 
y2 <- mean(babies$bwt[babies$smoke == 1], na.rm = TRUE)  
lines(c(0.8,1.2),rep(y1,2),lty=2,col=2, lwd=2)
lines(c(1.8,2.2),rep(y2,2),lty=2,col=2, lwd=2)


## ----q5b----------------------------------------------------------------------------------------------
# Filter the data to include only rows where smoke is 0 or 1
t_test <- t.test(bwt ~ smoke, data = babies)
print(t_test)


## ----q5c----------------------------------------------------------------------------------------------
set.seed(123)  
permute_t_stat <- function(data) {
  perm_smoke <- sample(data$smoke)  
  perm_t_test <- t.test(data$bwt ~ perm_smoke)  
  return(perm_t_test$statistic)  
}
perm_t_stat <- permute_t_stat(babies)
print(perm_t_stat)


## ----q5d----------------------------------------------------------------------------------------------
n_permutations <- 1000
perm_t_stats <- replicate(n_permutations, permute_t_stat(babies))
obs_t_stat <- t_test$statistic
p_value <- mean(abs(perm_t_stats) >= abs(obs_t_stat))
print(paste("P-value:", p_value))
plot(density(perm_t_stats), main = "Density of Permuted t-statistics", xlab = "t-statistics")
abline(v = obs_t_stat, col = "red", lwd = 2) 


## ----q5e----------------------------------------------------------------------------------------------
n_permutations <- 1000
perm_t_stats <- replicate(n_permutations, permute_t_stat(babies))
obs_t_stat <- t_test$statistic
p_value <- mean(abs(perm_t_stats) >= abs(obs_t_stat))
print(paste("P-value:", p_value))
plot(density(perm_t_stats), main = "Density of Permuted t-statistic", xlab = "Permuted t-statistic")
abline(v = obs_t_stat, col = "red", lwd = 2) 
curve(dt(x, df = t_test$parameter), add = TRUE, col = "blue", lwd = 2)
legend("topright", legend = c("Permuted t-stats", "Theoretical t-distribution"), 
       col = c("black", "blue"), lwd = 2)