<jupyter_start><jupyter_text>Confidence Interval SimulationThis notebook works through a simulation on confidence intervals, letting us see firsthand where the randomness is coming from.   PreambleFirst we load some useful packages and set the seed for our random number generator (this chooses a starting point, allowing us all the get the same results when starting with the same seed).<jupyter_code>set.seed(3482)
library(tidyverse)<jupyter_output><empty_output><jupyter_text>Simulating A Sample of size n=50Here we begin by taking a sample of size 50 for $X \sim N(4,4)$ and construct a confidence interval.<jupyter_code># One sample CI
n <- 50
# sample of 50 obs drawn from an N(4,4) distribution
sample <- rnorm(n,mean = 4,sd = 2)

# get the mean and sample standard deviation
m <- mean(sample)
s <- sd(sample)

# Calculate lower and upper bounds for a 95% confidence interval
ci_low <- m - 1.96*s/sqrt(n)
ci_upper <- m + 1.96*s/sqrt(n) 

one_sim <- data.frame(Mean = m,
                  Standard.Deviation = s,
                  CI_Lower = ci_low,
                  CI_Upper = ci_upper)
one_sim<jupyter_output><empty_output><jupyter_text>Simulating 1,000 samples of size n = 50Here we build a function to draw 1,000 samples of size 50 and construct a 95\% confidence interval for all samples.<jupyter_code># Function to perform multiple draws
simulate_samples <- function(nobs, nsim){
  res <- data.frame(index = 1:nsim,
                    m = 1:nsim,
                    s = 1:nsim)
  for(i in 1:nsim){
    sample <- rnorm(nobs,mean = 4,sd = 2)
    res$m[i] <- mean(sample)
    res$s[i] <- sd(sample)
    
  }
  return(res)
}

samples <- simulate_samples(nobs = n, nsim = 1000)
head(samples)
tail(samples)<jupyter_output><empty_output><jupyter_text>We now have 1,000 different samples, each with their own sample mean and sample standart deviation. We can now construct a confidence interval for all these samples.<jupyter_code>samples <- samples %>%
    mutate(ci_low = m - 1.96*s/sqrt(n),
           ci_upper = m + 1.96*s/sqrt(n))
head(samples)
tail(samples)<jupyter_output><empty_output><jupyter_text>Failures of our Confidence IntervalsNow let's check to see how many of these confidence intervals __do not__ include our true population mean, $\mu = 4$. We'll do this by creating binary variables that equal to one if the lower bound is greater than 4, or if the upper bound is less than 4. To get a sense of how often the confidence interval doesn't contain the true mean of 4, we can sum up the positive values for our $bad_lower$ and $bad_upper$ variables and divide by the number of simulations.<jupyter_code>samples <- samples %>%
    mutate(bad_lower = (4 < ci_low),
           bad_upper = (ci_upper < 4))
head(samples)
tail(samples)

bad_frac <- (sum(samples$bad_lower) + sum(samples$bad_upper))/ 1000
bad_frac<jupyter_output><empty_output><jupyter_text>And we see that the 95% confidence intervals for our simulations do not contain the true population mean in 5.2\% of the samples, almost exactly equal to 5%!   Plotting the Distribution of $\bar{X}$The below plot shows our distribution for $\bar{X}$. The red lines represent the upper and lower bounds for a 95% confidence interval using the true _population_ standard deviation. Since we generated our samples from an $N(4,4)$ distribution, we know that $\sigma = 2$ and can calculate the confidence interval as$$ \text{95% True Confidence Interval} ~=~ \left[\bar{x} - 1.96 \cdot \frac{\sigma}{\sqrt{n}},~ \bar{x} + 1.96 \cdot \frac{\sigma}{\sqrt{n}}\right]$$$$ \text{95% True Confidence Interval} ~=~ \left[\bar{x} - 1.96 \cdot \frac{2}{\sqrt{50}},~ \bar{x} + 1.96 \cdot \frac{2}{\sqrt{50}}\right]$$$$ \text{95% True Confidence Interval} ~=~ \left[3.4456, 4.5544 \right]$$<jupyter_code>true_low <- 4-1.96*2/sqrt(n)
true_upper <- 4+1.96*2/sqrt(n)

ggplot(samples, aes(m)) +
  geom_density() +
  geom_vline(xintercept = true_low, colour = "red", alpha = 0.6) +
  geom_vline(xintercept = true_upper, colour = "red", alpha = 0.6) + 
  labs( title = "Distribution of Simulated Sample Means (1,000 draws of n=50)",
        x = "Mean",
        y = "Density")<jupyter_output><empty_output><jupyter_text>And we see three things:1. Our distribution looks very similar to a normal distribution. If we took more draws or increased our sample size, it would likely smooth out even more.2. The distribution is centered around the true mean $\mu = 4$3. It seems like a reasonable guess that the lower and upper confidence interval lines look like roughly 0.025 of the distribution's area fall beyond them  Plotting the Simulated Confidence IntervalsAnother thing we expect due to the symmetry of the normal distribution is that our confidence interval failures should be symmetric. That is, $bad\_lower$ and $bad\_upper$ should be about equal in number.Another way to evaluate this is by plotting the lower versus upper limits and comparing their shape to the horizontal and vertical lines at $\mu =4$. We would like to see a scatterplot that is symmetric around the 45 degree line through the point $(4,4)$.<jupyter_code>ggplot(samples) +
geom_point(aes(x = ci_low, y = ci_upper), shape = 1, alpha = 0.7) +
geom_abline(intercept = 8, slope = -1, linetype = "dotted") + 
geom_hline(yintercept = 4, colour = "red", alpha = 0.6) +
geom_vline(xintercept = 4, colour = "red", alpha = 0.6) +
labs(x = "Lower End of Confidence Interval",
    y = "Upper End of Confidence Interval",
    title = "Symmetry of Confidence Interval Bounds (1,000 draws of n=50)")<jupyter_output><empty_output><jupyter_text>What Happens When We Increase the Sample Size?Previously all our samples were of size $n=50$. What if we instead took samples of size $n' = 500$?<jupyter_code>nprime <- 500
samples_prime <- simulate_samples(nobs = nprime, nsim = 1000) %>%
    mutate(ci_low = m - 1.96*s/sqrt(nprime),
           ci_upper = m + 1.96*s/sqrt(nprime))

head(samples_prime)
tail(samples_prime)<jupyter_output><empty_output><jupyter_text>First off, we can see that the sample means and standard deviations are much more tightly clustered around the true population mean (4) and standard deviation (2). As a result, our 95% confidence intervals are also much tighter around the true mean.<jupyter_code>samples_prime <- samples_prime %>%
    mutate(bad_lower = (4 < ci_low),
           bad_upper = (ci_upper < 4))
head(samples_prime)
tail(samples_prime)

bad_frac_prime <- (sum(samples_prime$bad_lower) + sum(samples_prime$bad_upper))/ 1000
bad_frac_prime<jupyter_output><empty_output><jupyter_text>And once again we have about 5% of the (now much narrower) confidence intervals not including the true mean $\mu =4$.<jupyter_code>true_low_prime <- 4-1.96*2/sqrt(nprime)
true_upper_prime <- 4+1.96*2/sqrt(nprime)

ggplot(samples_prime, aes(m)) +
  geom_density() +
  geom_vline(xintercept = true_low_prime, colour = "red", alpha = 0.6) +
  geom_vline(xintercept = true_upper_prime, colour = "red", alpha = 0.6) + 
  labs( title = "Distribution of Simulated Sample Means (1,000 draws of n=500)",
        x = "Mean",
        y = "Density")<jupyter_output><empty_output><jupyter_text>Our distribution now looks a lot more like a normal distribution, looking much smoother overall. Notice how it is also much more tightly concentrated around the true population mean.<jupyter_code>ggplot(samples_prime) +
geom_point(aes(x = ci_low, y = ci_upper), shape = 1, alpha = 0.7) +
geom_abline(intercept = 8, slope = -1, linetype = "dotted") + 
geom_hline(yintercept = 4, colour = "red", alpha = 0.6) +
geom_vline(xintercept = 4, colour = "red", alpha = 0.6) +
labs(x = "Lower End of Confidence Interval",
    y = "Upper End of Confidence Interval",
    title = "Symmetry of Confidence Interval Bounds (1,000 draws of n=500)")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10778135048231512, 'fraction_numerical': 0.029967845659163987, 'mean_word_length': 4.39251040221914, 'pattern_counts': {'":': 0, '<': 73, '<?xml version=': 0, '>': 49, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '428552', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2672, 'n_tokens_neox': 2440, 'n_words': 890, 'path': 'Summer-19/Simulations/Confidence Interval Simulation.ipynb', 'repo_name': 'sberto2019/ENVECON-118', 'size': '7775'}