library(Stat2Data)
library(ggplot2)
data(Diamonds)

# Step 1: Population Parameters
mu <- mean(Diamonds$TotalPrice)
sigma <- sd(Diamonds$TotalPrice)

# Step 2: Simulation (n = 30, 1000 reps)
n <- 30
reps <- 1000
set.seed(123) # Ensures reproducible random sampling

sample_means <- replicate(reps, {
  sample_data <- sample(Diamonds$TotalPrice, size = n, replace = TRUE)
  mean(sample_data)
})

sim_data <- data.frame(sample_means = sample_means)

# Step 3: Theoretical Standard Error & Visualization
SE <- sigma / sqrt(n)

ggplot(sim_data, aes(x = sample_means)) +
  geom_histogram(aes(y = after_stat(density)), color = "white", fill = "salmon", bins = 30) +
  stat_function(fun = dnorm, args = list(mean = mu, sd = SE), color = "blue", size = 1) +
  labs(title = "Sampling Distribution of Mean Diamond Prices (n = 30)",
       x = "Sample Mean Total Price ($)",
       y = "Density") +
  theme_minimal()

# Step 4: Probability Calculation (P(X_bar > 10000))
prob_theoretical <- 1 - pnorm(10000, mean = mu, sd = SE)
prob_empirical <- mean(sample_means > 10000)

cat("Theoretical P(Mean > 10000):", prob_theoretical, "\n")

## Theoretical P(Mean > 10000): 0.03632525

cat("Empirical P(Mean > 10000):  ", prob_empirical, "\n\n")

## Empirical P(Mean > 10000):   0.05

# Step 5: Cutoff Calculation (90% above / 10th percentile)
cutoff_theoretical <- qnorm(0.10, mean = mu, sd = SE)
cutoff_empirical <- quantile(sample_means, 0.10)

cat("Theoretical 10th percentile cutoff:", cutoff_theoretical, "\n")

## Theoretical 10th percentile cutoff: 5629.452

cat("Empirical 10th percentile cutoff:  ", cutoff_empirical, "\n")

## Empirical 10th percentile cutoff:   5816.595