# ====================================================================
# Week 6: Probability and Statistical Foundations
# Kmex Consult - R Data Science Course
# Open in RStudio and run each block with Ctrl+Enter (Cmd+Enter on Mac).
# ====================================================================

# Understanding randomness, distributions, and inference

# ====================================================================
# WEEK 6 OVERVIEW
# ====================================================================
# This week is FOUNDATIONAL - everything from Weeks 7-12 builds on these concepts.
# WHAT YOU NEED:
# - Understanding of probability
# - Knowledge of distributions
# - Concept of sampling
# - Central Limit Theorem
# - Confidence intervals
# - Basics of hypothesis testing (detailed in Week 7)
# KEY CONCEPTS:
# Probability: Chance something happens (0 to 1)
# Distribution: Pattern of values
# Population: All data (usually unknown)
# Sample: Subset we can observe
# Inference: Drawing conclusions about population from sample
# WHY THIS MATTERS:
# - All statistical tests depend on these
# - Understanding uncertainty is critical
# - Makes you different from analysts who just run tests
# - Enables you to choose right test
# - Helps interpret results correctly
# THIS WEEK YOU'LL LEARN:
# ✓ Probability fundamentals and rules
# ✓ Bayes' theorem (probability updating)
# ✓ Discrete distributions (binomial, Poisson)
# ✓ Continuous distributions (normal, uniform, exponential)
# ✓ Sampling and sampling distributions
# ✓ Central Limit Theorem
# ✓ Confidence intervals
# ✓ Sampling procedures for data collection

# ====================================================================
# SESSION 1: PROBABILITY AND DISTRIBUTIONS
# ====================================================================

# --- Duration: 2 hours ---

# --- 1.1 Probability Fundamentals ---
# PROBABILITY - FUNDAMENTAL CONCEPT:
# Definition: Probability of event A = Number of favorable outcomes / Total outcomes


# Coin flip
P(heads) = 1 / 2 = 0.5

# Die roll
P(rolling 3) = 1 / 6 = 0.167

# Card from deck
P(ace) = 4 / 52 = 0.077

# PROBABILITY RULES:


# Rule 1: Probability between 0 and 1
P(A) >= 0 and P(A) <= 1

# Rule 2: All outcomes sum to 1
P(A) + P(not A) = 1

# Rule 3: Addition rule (A OR B)
P(A or B) = P(A) + P(B) - P(A and B)

# Rule 4: Multiplication rule (A AND B)
P(A and B) = P(A) * P(B | A) # If independent: P(A) * P(B)

# Rule 5: Conditional probability
P(A | B) = P(A and B) / P(B) # Probability of A given B happened

# EXAMPLE: Disease testing
P(disease) = 0.01
P(positive | disease) = 0.99 # Sensitivity
P(positive | no disease) = 0.05 # False positive rate

# What's P(disease | positive)?
# P(disease | positive) = P(positive | disease) * P(disease) / P(positive)
# Requires Bayes' theorem...

# BAYES' THEOREM - UPDATING BELIEFS:


# P(A | B) = P(B | A) * P(A) / P(B)

# Or: Updated belief = Likelihood * Prior / Evidence

# EXAMPLE: Email spam detection
# Prior: 20% of emails are spam
# Likelihood: 90% of spam have "viagra"
# Prior: 10% of legitimate have "viagra"
# Question: Email has "viagra", probability it's spam?

prior_spam = 0.20
likelihood_given_spam = 0.90
likelihood_given_not_spam = 0.10
prior_not_spam = 0.80

# Calculate evidence (total probability of "viagra")
evidence = (likelihood_given_spam * prior_spam) +
(likelihood_given_not_spam * prior_not_spam)
# = (0.90 * 0.20) + (0.10 * 0.80) = 0.18 + 0.08 = 0.26

# Apply Bayes
posterior = (likelihood_given_spam * prior_spam) / evidence
# = (0.90 * 0.20) / 0.26 = 0.18 / 0.26 = 0.69 (69% spam)

# In R:
bayes_theorem <- function(prior, likelihood_yes, likelihood_no) {
evidence <- (likelihood_yes * prior) + (likelihood_no * (1 - prior))
posterior <- (likelihood_yes * prior) / evidence
posterior
}

spam_prob <- bayes_theorem(0.20, 0.90, 0.10)
print(spam_prob) # [1] 0.6923077

# DISCRETE DISTRIBUTIONS:
# BINOMIAL DISTRIBUTION (Success/Failure, n trials):


# Parameters: n (number trials), p (probability success)

# Example: 10 coin flips, P(heads) = 0.5
dbinom(5, size=10, prob=0.5) # Probability of exactly 5 heads
# [1] 0.2460938

pbinom(5, size=10, prob=0.5) # Probability of ≤ 5 heads (cumulative)
# [1] 0.6230469

rbinom(10, size=10, prob=0.5) # Simulate 10 coin flips

# Plot distribution
library(ggplot2)
x <- 0:10
y <- dbinom(x, size=10, prob=0.5)
ggplot(tibble(x, y), aes(x, y)) +
geom_bar(stat="identity", fill="steelblue") +
labs(title="Binomial Distribution (n=10, p=0.5)",
x="Number of Successes", y="Probability")

# POISSON DISTRIBUTION (Rare events over time):


# Parameter: lambda (mean rate)

# Example: Customer service calls per hour (average 3)
dpois(2, lambda=3) # Probability of exactly 2 calls
# [1] 0.224042

ppois(5, lambda=3) # Probability of ≤ 5 calls
# [1] 0.9161157

rpois(10, lambda=3) # Simulate 10 hours of calls

# CONTINUOUS DISTRIBUTIONS:
# NORMAL DISTRIBUTION (Bell curve - most important!):


# Parameters: mean (μ), sd (σ)

# Standard normal (mean=0, sd=1)
dnorm(0) # Density at 0
# [1] 0.3989423

pnorm(0) # Probability ≤ 0
# [1] 0.5

pnorm(1.96) # Probability ≤ 1.96
# [1] 0.975 (95th percentile)

qnorm(0.975) # Value at 97.5th percentile
# [1] 1.959964

rnorm(1000) # Generate 1000 normal values

# Visualize
x <- seq(-4, 4, by=0.1)
y <- dnorm(x)
ggplot(tibble(x, y), aes(x, y)) +
geom_line(color="steelblue", size=1) +
labs(title="Normal Distribution (μ=0, σ=1)",
x="Value", y="Density")

# 68-95-99.7 rule
pnorm(1) - pnorm(-1) # ~68% within 1 SD
# [1] 0.6826895

pnorm(2) - pnorm(-2) # ~95% within 2 SD
# [1] 0.9544997

pnorm(3) - pnorm(-3) # ~99.7% within 3 SD
# [1] 0.9973002

# UNIFORM DISTRIBUTION (All values equally likely):


# Rectangle shape

dunif(0.5, min=0, max=1) # Density
punif(0.5, min=0, max=1) # CDF
runif(1000, min=0, max=1) # Generate values

# EXPONENTIAL DISTRIBUTION (Time to event):


# Parameter: rate (λ)

dexp(1, rate=1)
pexp(1, rate=1)
rexp(1000, rate=1)


# --- 1.2 Sampling and Central Limit Theorem ---
# SAMPLING DISTRIBUTION - KEY CONCEPT:
# Imagine repeating your experiment 1000 times:
# - Each time you get a different sample mean
# - These means form a distribution (sampling distribution!)
# - This distribution has its own mean and standard error


# Simulate
set.seed(42)
sample_means <- vector("numeric", 1000)

for (i in 1:1000) {
sample <- rnorm(30, mean=100, sd=15)
sample_means[i] <- mean(sample)
}

# Distribution of sample means
mean(sample_means) # ~100 (same as population)
sd(sample_means) # ~2.74 (much smaller!)

# Plot
ggplot(tibble(sample_means), aes(x=sample_means)) +
geom_histogram(bins=30, fill="steelblue") +
labs(title="Sampling Distribution of Mean",
x="Sample Mean", y="Frequency")

# CENTRAL LIMIT THEOREM - MOST IMPORTANT THEOREM:
# "If we repeatedly sample from ANY population, the distribution of sample means is approximately normal"
# IMPLICATIONS:
# 1. Works regardless of original distribution shape
# 2. Sample mean ≈ population mean
# 3. Standard error = population SD / sqrt(n)
# 4. Larger n → smaller standard error (more precise)


# Demonstrate with non-normal data (exponential)
set.seed(42)

# Original is NOT normal
original <- rexp(10000, rate=1)
ggplot(tibble(original), aes(x=original)) +
geom_histogram(bins=50, fill="coral") +
labs(title="Original Distribution (Exponential - NOT Normal)")

# BUT sample means ARE normal!
sample_means <- vector("numeric", 1000)
for (i in 1:1000) {
sample <- rexp(30, rate=1)
sample_means[i] <- mean(sample)
}

ggplot(tibble(sample_means), aes(x=sample_means)) +
geom_histogram(bins=30, fill="steelblue") +
geom_density(aes(y=..density..), color="red") +
labs(title="Distribution of Sample Means (NORMAL!)")

# Shapiro-Wilk test for normality
shapiro.test(sample_means) # p > 0.05 → normal

# STANDARD ERROR:
# Standard Error = SD / sqrt(n)


sd_population <- 15
n_sample <- 30
se <- sd_population / sqrt(n_sample)
print(se) # [1] 2.738613

# Larger samples → smaller standard error
n_values <- c(10, 30, 100, 500)
se_values <- 15 / sqrt(n_values)

ggplot(tibble(n = n_values, se = se_values), aes(n, se)) +
geom_line() +
geom_point() +
labs(title="Standard Error vs Sample Size",
x="Sample Size", y="Standard Error")

# CONFIDENCE INTERVALS - UNCERTAINTY QUANTIFICATION:
# What's a 95% confidence interval?
# "If we repeated our experiment 100 times, 95% of CIs would contain true parameter"
# NOT: "95% probability true value is in this interval" (common misconception)


# Calculate 95% CI
sample <- rnorm(100, mean=100, sd=15)
sample_mean <- mean(sample)
se <- sd(sample) / sqrt(100)
margin_error <- 1.96 * se # 1.96 for 95%

ci_lower <- sample_mean - margin_error
ci_upper <- sample_mean + margin_error

print(paste(ci_lower, "to", ci_upper))

# Using t-distribution (better when n < 30)
t_critical <- qt(0.975, df=99) # 0.975 for 95% two-tailed
margin_error_t <- t_critical * se
ci_lower_t <- sample_mean - margin_error_t
ci_upper_t <- sample_mean + margin_error_t

# Using built-in functions
t.test(sample)$conf.int # Automatic!
# [1] 97.34 102.97

# Interpretation: We're 95% confident true mean is between 97.34 and 102.97
# NOT: 95% probability true mean is in this range (wrong!)

# FACTORS AFFECTING CONFIDENCE INTERVALS:


# Larger confidence level → wider interval
# 90% CI narrower than 95% CI narrower than 99% CI

# Larger sample → narrower interval
# More data = more precise estimate

# Larger SD → wider interval
# More variability = less precise

# Visualization
df_example <- tibble(
ci_level = c("90%", "95%", "99%"),
lower = c(-1.645, -1.960, -2.576),
upper = c(1.645, 1.960, 2.576)
)

ggplot(df_example, aes(y=ci_level)) +
geom_errorbarh(aes(xmin=lower, xmax=upper), height=0.2) +
geom_point(aes(x=0)) +
labs(title="Confidence Intervals at Different Levels",
x="Value", y="Confidence Level")