{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "29b2c039",
   "metadata": {},
   "source": [
    "# Week 6: Probability, Distributions and Statistical Inference\n",
    "\n",
    "*Kmex Consult — Data Science Fundamentals Course*\n",
    "\n",
    "Work through each code cell in order. Edit and re-run to experiment."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d92c1b6",
   "metadata": {},
   "source": [
    "*Data Science Fundamentals Course*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5aae2fc2",
   "metadata": {},
   "source": [
    "## Week 6 Overview\n",
    "Probability and statistics form the theoretical foundation of data science. Understanding these concepts is essential for:\n",
    "- Interpreting data correctly\n",
    "- Building and evaluating machine learning models\n",
    "- Making data-driven decisions\n",
    "- Quantifying uncertainty\n",
    "- Conducting hypothesis tests\n",
    "\n",
    "This week focuses on fundamental concepts that underpin all statistical analysis and machine learning. You'll learn how to work with probability distributions, understand the Central Limit Theorem, calculate confidence intervals, and use simulations to understand complex statistical concepts.\n",
    "\n",
    "By the end of Week 6, you will be able to:\n",
    "- Understand and calculate probabilities\n",
    "- Work with discrete and continuous probability distributions\n",
    "- Apply the Central Limit Theorem\n",
    "- Calculate and interpret confidence intervals\n",
    "- Understand sampling distributions\n",
    "- Use bootstrapping for statistical estimation\n",
    "- Understand the difference between Bayesian and Frequentist approaches\n",
    "- Use simulations to understand statistical concepts\n",
    "\n",
    "Week 6 is divided into three 2-hour sessions:\n",
    "- Session 1: Probability Fundamentals and Distributions\n",
    "- Session 2: Sampling Distributions and the Central Limit Theorem\n",
    "- Session 3: Confidence Intervals, Bootstrapping, and Bayesian Thinking"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4106f21",
   "metadata": {},
   "source": [
    "## SESSION 1: Probability Fundamentals and Distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "003611de",
   "metadata": {},
   "source": [
    "### Duration: 2 hours"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f29d7e0c",
   "metadata": {},
   "source": [
    "### 1.1 Probability Fundamentals\n",
    "Probability is the mathematical study of randomness and uncertainty. It quantifies how likely an event is to occur.\n",
    "\n",
    "Key concepts:\n",
    "- Sample Space: All possible outcomes of an experiment\n",
    "- Event: A subset of the sample space\n",
    "- Probability: A number between 0 and 1 indicating likelihood\n",
    "- Mutually Exclusive: Events that cannot occur together\n",
    "- Independent: Events where one doesn't affect the other"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7879ab5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Basic probability\n",
    "# Probability = Number of favorable outcomes / Total possible outcomes\n",
    "\n",
    "print(\"=\"*50)\n",
    "print(\"BASIC PROBABILITY EXAMPLES\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "# Example 1: Rolling a die\n",
    "# Probability of rolling a 4\n",
    "favorable_outcomes = 1\n",
    "total_outcomes = 6\n",
    "prob_four = favorable_outcomes / total_outcomes\n",
    "print(f\"\\nP(rolling a 4) = {favorable_outcomes}/{total_outcomes} = {prob_four:.3f} or {prob_four*100:.1f}%\")\n",
    "\n",
    "# Example 2: Drawing from a deck\n",
    "# Probability of drawing a red card\n",
    "red_cards = 26\n",
    "total_cards = 52\n",
    "prob_red = red_cards / total_cards\n",
    "print(f\"P(red card) = {red_cards}/{total_cards} = {prob_red:.3f} or {prob_red*100:.1f}%\")\n",
    "\n",
    "# Example 3: Probability rules\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"PROBABILITY RULES\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "# Addition rule: P(A or B) = P(A) + P(B) - P(A and B)\n",
    "P_A = 0.3 # Probability of event A\n",
    "P_B = 0.4 # Probability of event B\n",
    "P_A_and_B = 0.1 # Probability of both occurring\n",
    "\n",
    "P_A_or_B = P_A + P_B - P_A_and_B\n",
    "print(f\"\\nAddition Rule: P(A or B) = P(A) + P(B) - P(A and B)\")\n",
    "print(f\"P(A or B) = {P_A} + {P_B} - {P_A_and_B} = {P_A_or_B}\")\n",
    "\n",
    "# Multiplication rule: P(A and B) = P(A) * P(B|A)\n",
    "# For independent events: P(A and B) = P(A) * P(B)\n",
    "P_heads = 0.5\n",
    "P_tails = 0.5\n",
    "P_two_heads = P_heads * P_heads\n",
    "print(f\"\\nMultiplication Rule (independent): P(two heads) = 0.5 × 0.5 = {P_two_heads}\")\n",
    "\n",
    "# Complement rule: P(not A) = 1 - P(A)\n",
    "P_rain = 0.3\n",
    "P_no_rain = 1 - P_rain\n",
    "print(f\"\\nComplement Rule: P(no rain) = 1 - P(rain) = 1 - {P_rain} = {P_no_rain}\")\n",
    "\n",
    "# Conditional probability: P(A|B) = P(A and B) / P(B)\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"CONDITIONAL PROBABILITY\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "# Example: Disease test accuracy\n",
    "# P(Disease|Positive Test)\n",
    "P_positive_given_disease = 0.95 # Sensitivity\n",
    "P_negative_given_no_disease = 0.90 # Specificity\n",
    "P_disease = 0.01 # Prior probability of disease\n",
    "\n",
    "# Using Bayes' theorem\n",
    "P_positive = (P_positive_given_disease * P_disease) + ((1 - P_negative_given_no_disease) * (1 - P_disease))\n",
    "P_disease_given_positive = (P_positive_given_disease * P_disease) / P_positive\n",
    "\n",
    "print(f\"Sensitivity (P(+|disease)): {P_positive_given_disease}\")\n",
    "print(f\"Specificity (P(-|no disease)): {P_negative_given_no_disease}\")\n",
    "print(f\"P(disease) prior: {P_disease}\")\n",
    "print(f\"P(+) overall: {P_positive:.4f}\")\n",
    "print(f\"P(disease|+) posterior: {P_disease_given_positive:.4f}\")\n",
    "print(f\"\\nInterpretation: Even with positive test, only {P_disease_given_positive*100:.1f}% probability of having disease!\")\n",
    "\n",
    "# Simulation: Rolling dice\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"SIMULATION: ROLLING DICE\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "np.random.seed(42)\n",
    "num_rolls = 100000\n",
    "\n",
    "# Simulate rolling a die\n",
    "rolls = np.random.randint(1, 7, num_rolls)\n",
    "\n",
    "# Calculate probabilities\n",
    "for face in range(1, 7):\n",
    "prob = (rolls == face).sum() / num_rolls\n",
    "print(f\"P(rolling {face}) = {prob:.4f}\")\n",
    "\n",
    "print(f\"\\nExpected: {1/6:.4f} for each face\")\n",
    "\n",
    "# Simulate coin flips\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"SIMULATION: COIN FLIPS\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "flips = np.random.randint(0, 2, 100000)\n",
    "prob_heads = (flips == 1).sum() / 100000\n",
    "prob_tails = (flips == 0).sum() / 100000\n",
    "\n",
    "print(f\"P(Heads) = {prob_heads:.4f}\")\n",
    "print(f\"P(Tails) = {prob_tails:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a7a4849",
   "metadata": {},
   "source": [
    "### 1.2 Discrete Probability Distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31e43fe2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print(\"=\"*50)\n",
    "print(\"DISCRETE DISTRIBUTIONS\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "# 1. BERNOULLI DISTRIBUTION\n",
    "# Single trial with two outcomes (success/failure)\n",
    "print(\"\\n1. BERNOULLI DISTRIBUTION\")\n",
    "print(\"Example: Single coin flip, single product working/failing\")\n",
    "\n",
    "p = 0.7 # Probability of success\n",
    "bernoulli = stats.bernoulli(p)\n",
    "\n",
    "print(f\"P(success) = {p}\")\n",
    "print(f\"P(failure) = {1-p}\")\n",
    "print(f\"Mean = {bernoulli.mean()}\")\n",
    "print(f\"Variance = {bernoulli.var()}\")\n",
    "\n",
    "# 2. BINOMIAL DISTRIBUTION\n",
    "# Number of successes in n independent trials\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"2. BINOMIAL DISTRIBUTION\")\n",
    "print(\"Example: Number of heads in 10 coin flips, defects in 100 items\")\n",
    "\n",
    "n = 10 # Number of trials\n",
    "p = 0.5 # Probability of success per trial\n",
    "binomial = stats.binom(n, p)\n",
    "\n",
    "print(f\"\\nParameters: n={n} trials, p={p}\")\n",
    "print(f\"Mean = {binomial.mean()}\")\n",
    "print(f\"Variance = {binomial.var()}\")\n",
    "\n",
    "# Probability of exactly 5 successes\n",
    "prob_5 = binomial.pmf(5)\n",
    "print(f\"P(X=5) = {prob_5:.4f}\")\n",
    "\n",
    "# Probability of 5 or fewer successes\n",
    "prob_le5 = binomial.cdf(5)\n",
    "print(f\"P(X≤5) = {prob_le5:.4f}\")\n",
    "\n",
    "# Simulate\n",
    "simulated = np.random.binomial(n, p, 10000)\n",
    "print(f\"\\nSimulated mean: {simulated.mean():.3f}\")\n",
    "print(f\"Actual mean: {binomial.mean():.3f}\")\n",
    "\n",
    "# 3. POISSON DISTRIBUTION\n",
    "# Number of events in fixed interval\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"3. POISSON DISTRIBUTION\")\n",
    "print(\"Example: Number of customer arrivals per hour, defects per batch\")\n",
    "\n",
    "lambda_param = 3 # Average number of events\n",
    "poisson = stats.poisson(lambda_param)\n",
    "\n",
    "print(f\"\\nParameter (λ) = {lambda_param}\")\n",
    "print(f\"Mean = {poisson.mean()}\")\n",
    "print(f\"Variance = {poisson.var()}\")\n",
    "\n",
    "# Probability of exactly 5 events\n",
    "prob_5 = poisson.pmf(5)\n",
    "print(f\"P(X=5) = {prob_5:.4f}\")\n",
    "\n",
    "# Probability of 5 or fewer events\n",
    "prob_le5 = poisson.cdf(5)\n",
    "print(f\"P(X≤5) = {prob_le5:.4f}\")\n",
    "\n",
    "# Simulate\n",
    "simulated = np.random.poisson(lambda_param, 10000)\n",
    "print(f\"\\nSimulated mean: {simulated.mean():.3f}\")\n",
    "print(f\"Actual mean: {poisson.mean():.3f}\")\n",
    "\n",
    "# Visualization\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
    "\n",
    "# Binomial\n",
    "x = range(0, n+1)\n",
    "pmf = binomial.pmf(x)\n",
    "axes[0].bar(x, pmf, alpha=0.7, color='blue')\n",
    "axes[0].set_title(f'Binomial(n={n}, p={p})')\n",
    "axes[0].set_xlabel('Number of successes')\n",
    "axes[0].set_ylabel('Probability')\n",
    "axes[0].grid(True, alpha=0.3)\n",
    "\n",
    "# Poisson\n",
    "x = range(0, 10)\n",
    "pmf = poisson.pmf(x)\n",
    "axes[1].bar(x, pmf, alpha=0.7, color='green')\n",
    "axes[1].set_title(f'Poisson(λ={lambda_param})')\n",
    "axes[1].set_xlabel('Number of events')\n",
    "axes[1].set_ylabel('Probability')\n",
    "axes[1].grid(True, alpha=0.3)\n",
    "\n",
    "# Distribution comparison\n",
    "axes[2].bar([0, 1], [1-p, p], alpha=0.7, color='red', label='Bernoulli')\n",
    "axes[2].set_title(f'Bernoulli(p={p})')\n",
    "axes[2].set_xlabel('Outcome')\n",
    "axes[2].set_ylabel('Probability')\n",
    "axes[2].set_xticks([0, 1])\n",
    "axes[2].set_xticklabels(['Failure', 'Success'])\n",
    "axes[2].grid(True, alpha=0.3)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "# Real-world example: Product quality\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"REAL-WORLD: PRODUCT QUALITY CONTROL\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "# Defect rate is 2%, inspecting 50 units\n",
    "n_units = 50\n",
    "defect_rate = 0.02\n",
    "quality = stats.binom(n_units, defect_rate)\n",
    "\n",
    "print(f\"\\nDefect rate: {defect_rate*100}%\")\n",
    "print(f\"Inspecting: {n_units} units\")\n",
    "print(f\"Expected defects: {quality.mean()}\")\n",
    "print(f\"P(0 defects) = {quality.pmf(0):.4f}\")\n",
    "print(f\"P(≥2 defects) = {1 - quality.cdf(1):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4563191",
   "metadata": {},
   "source": [
    "### 1.3 Continuous Probability Distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a268cb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print(\"=\"*50)\n",
    "print(\"CONTINUOUS DISTRIBUTIONS\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "# 1. NORMAL (GAUSSIAN) DISTRIBUTION\n",
    "print(\"\\n1. NORMAL DISTRIBUTION\")\n",
    "print(\"Example: Heights, test scores, measurement errors\")\n",
    "\n",
    "mu = 100 # Mean\n",
    "sigma = 15 # Standard deviation\n",
    "normal = stats.norm(mu, sigma)\n",
    "\n",
    "print(f\"\\nParameters: μ={mu}, σ={sigma}\")\n",
    "print(f\"Mean = {normal.mean()}\")\n",
    "print(f\"Variance = {normal.var()}\")\n",
    "\n",
    "# Probabilities\n",
    "prob_less_100 = normal.cdf(100)\n",
    "print(f\"P(X < 100) = {prob_less_100:.4f}\")\n",
    "\n",
    "prob_100_110 = normal.cdf(110) - normal.cdf(100)\n",
    "print(f\"P(100 < X < 110) = {prob_100_110:.4f}\")\n",
    "\n",
    "# Percentiles\n",
    "percentile_90 = normal.ppf(0.90)\n",
    "print(f\"90th percentile = {percentile_90:.2f}\")\n",
    "\n",
    "# Z-scores (standardized)\n",
    "z_score = (110 - mu) / sigma\n",
    "print(f\"Z-score for X=110: {z_score:.3f}\")\n",
    "\n",
    "# Standard normal (μ=0, σ=1)\n",
    "print(f\"\\nP(Z < 0) in standard normal = {stats.norm.cdf(0):.4f}\")\n",
    "print(f\"P(-1 < Z < 1) = {stats.norm.cdf(1) - stats.norm.cdf(-1):.4f}\")\n",
    "print(f\"P(-2 < Z < 2) = {stats.norm.cdf(2) - stats.norm.cdf(-2):.4f}\")\n",
    "print(f\"P(-3 < Z < 3) = {stats.norm.cdf(3) - stats.norm.cdf(-3):.4f}\")\n",
    "\n",
    "# 2. EXPONENTIAL DISTRIBUTION\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"2. EXPONENTIAL DISTRIBUTION\")\n",
    "print(\"Example: Time until next event, equipment failure time\")\n",
    "\n",
    "lambda_param = 0.5 # Rate parameter\n",
    "exponential = stats.expon(scale=1/lambda_param)\n",
    "\n",
    "print(f\"\\nParameter (λ) = {lambda_param}\")\n",
    "print(f\"Mean = {exponential.mean():.3f}\")\n",
    "print(f\"Variance = {exponential.var():.3f}\")\n",
    "\n",
    "# Probability density\n",
    "x = np.linspace(0, 10, 100)\n",
    "pdf = exponential.pdf(x)\n",
    "\n",
    "# Cumulative probability\n",
    "prob_less_2 = exponential.cdf(2)\n",
    "print(f\"P(X < 2) = {prob_less_2:.4f}\")\n",
    "\n",
    "# 3. UNIFORM DISTRIBUTION\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"3. UNIFORM DISTRIBUTION\")\n",
    "print(\"Example: Random number generation, continuous random selection\")\n",
    "\n",
    "a, b = 0, 10 # Bounds\n",
    "uniform = stats.uniform(a, b-a)\n",
    "\n",
    "print(f\"\\nBounds: [{a}, {b}]\")\n",
    "print(f\"Mean = {uniform.mean():.3f}\")\n",
    "print(f\"Variance = {uniform.var():.3f}\")\n",
    "\n",
    "prob_3_7 = uniform.cdf(7) - uniform.cdf(3)\n",
    "print(f\"P(3 < X < 7) = {prob_3_7:.4f}\")\n",
    "\n",
    "# Visualization\n",
    "fig, axes = plt.subplots(2, 2, figsize=(12, 10))\n",
    "\n",
    "# Normal distribution\n",
    "x = np.linspace(mu - 4*sigma, mu + 4*sigma, 100)\n",
    "pdf = normal.pdf(x)\n",
    "axes[0, 0].plot(x, pdf, 'b-', linewidth=2, label='Normal')\n",
    "axes[0, 0].fill_between(x, pdf, alpha=0.3)\n",
    "axes[0, 0].axvline(mu, color='r', linestyle='--', label=f'Mean = {mu}')\n",
    "axes[0, 0].set_title(f'Normal Distribution μ={mu}, σ={sigma}')\n",
    "axes[0, 0].set_xlabel('Value')\n",
    "axes[0, 0].set_ylabel('Probability Density')\n",
    "axes[0, 0].legend()\n",
    "axes[0, 0].grid(True, alpha=0.3)\n",
    "\n",
    "# Exponential distribution\n",
    "x = np.linspace(0, 10, 100)\n",
    "pdf = exponential.pdf(x)\n",
    "axes[0, 1].plot(x, pdf, 'g-', linewidth=2)\n",
    "axes[0, 1].fill_between(x, pdf, alpha=0.3, color='green')\n",
    "axes[0, 1].set_title(f'Exponential Distribution λ={lambda_param}')\n",
    "axes[0, 1].set_xlabel('Value')\n",
    "axes[0, 1].set_ylabel('Probability Density')\n",
    "axes[0, 1].grid(True, alpha=0.3)\n",
    "\n",
    "# Uniform distribution\n",
    "x = np.linspace(a-1, b+1, 100)\n",
    "pdf = uniform.pdf(x)\n",
    "axes[1, 0].plot(x, pdf, 'r-', linewidth=2)\n",
    "axes[1, 0].fill_between(x, pdf, alpha=0.3, color='red')\n",
    "axes[1, 0].set_title(f'Uniform Distribution [{a}, {b}]')\n",
    "axes[1, 0].set_xlabel('Value')\n",
    "axes[1, 0].set_ylabel('Probability Density')\n",
    "axes[1, 0].grid(True, alpha=0.3)\n",
    "\n",
    "# Multiple normal distributions for comparison\n",
    "x = np.linspace(-50, 250, 100)\n",
    "for mu_i in [100, 150, 200]:\n",
    "pdf = stats.norm.pdf(x, mu_i, 15)\n",
    "axes[1, 1].plot(x, pdf, linewidth=2, label=f'μ={mu_i}')\n",
    "axes[1, 1].set_title('Normal Distributions with Different Means')\n",
    "axes[1, 1].set_xlabel('Value')\n",
    "axes[1, 1].set_ylabel('Probability Density')\n",
    "axes[1, 1].legend()\n",
    "axes[1, 1].grid(True, alpha=0.3)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "# Real-world example: Test scores\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"REAL-WORLD: TEST SCORE ANALYSIS\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "mean_score = 75\n",
    "std_score = 8\n",
    "test_dist = stats.norm(mean_score, std_score)\n",
    "\n",
    "print(f\"\\nMean score: {mean_score}\")\n",
    "print(f\"Standard deviation: {std_score}\")\n",
    "print(f\"P(score < 70) = {test_dist.cdf(70):.4f} or {test_dist.cdf(70)*100:.2f}%\")\n",
    "print(f\"P(score > 90) = {1 - test_dist.cdf(90):.4f} or {(1-test_dist.cdf(90))*100:.2f}%\")\n",
    "print(f\"P(70 < score < 85) = {test_dist.cdf(85) - test_dist.cdf(70):.4f} or {(test_dist.cdf(85) - test_dist.cdf(70))*100:.2f}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09e95faf",
   "metadata": {},
   "source": [
    "## SESSION 2: Sampling Distributions and Central Limit Theorem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dea30e5",
   "metadata": {},
   "source": [
    "### Duration: 2 hours"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d407ccbc",
   "metadata": {},
   "source": [
    "### 2.1 Sampling Distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8da6bb92",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print(\"=\"*50)\n",
    "print(\"SAMPLING DISTRIBUTIONS\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "KEY CONCEPTS:\n",
    "- Population: All items of interest\n",
    "- Sample: Subset of population\n",
    "- Sampling distribution: Distribution of sample statistics\n",
    "- Standard Error: Standard deviation of sampling distribution\n",
    "\"\"\")\n",
    "\n",
    "# Create a population\n",
    "np.random.seed(42)\n",
    "population = np.random.normal(100, 15, 100000)\n",
    "\n",
    "print(f\"\\nPopulation statistics:\")\n",
    "print(f\"Mean: {population.mean():.2f}\")\n",
    "print(f\"Std Dev: {population.std():.2f}\")\n",
    "\n",
    "# Draw many samples and calculate means\n",
    "sample_size = 30\n",
    "num_samples = 10000\n",
    "sample_means = []\n",
    "\n",
    "for i in range(num_samples):\n",
    "sample = np.random.choice(population, sample_size, replace=False)\n",
    "sample_means.append(sample.mean())\n",
    "\n",
    "sample_means = np.array(sample_means)\n",
    "\n",
    "print(f\"\\n\" + \"=\"*50)\n",
    "print(\"SAMPLING DISTRIBUTION OF MEANS\")\n",
    "print(\"=\"*50)\n",
    "print(f\"Sample size: {sample_size}\")\n",
    "print(f\"Number of samples: {num_samples}\")\n",
    "print(f\"\\nMean of sample means: {sample_means.mean():.2f}\")\n",
    "print(f\"Std Dev of sample means (Standard Error): {sample_means.std():.2f}\")\n",
    "\n",
    "# Compare to theoretical standard error\n",
    "se_theoretical = population.std() / np.sqrt(sample_size)\n",
    "print(f\"Theoretical SE: {se_theoretical:.2f}\")\n",
    "\n",
    "# Effect of sample size\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"EFFECT OF SAMPLE SIZE\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "sample_sizes = [10, 30, 100, 300]\n",
    "results = {}\n",
    "\n",
    "for n in sample_sizes:\n",
    "sample_means_n = []\n",
    "for i in range(5000):\n",
    "sample = np.random.choice(population, n, replace=False)\n",
    "sample_means_n.append(sample.mean())\n",
    "\n",
    "results[n] = np.array(sample_means_n)\n",
    "se = population.std() / np.sqrt(n)\n",
    "print(f\"Sample size {n:3d}: SE = {results[n].std():.2f} (theoretical: {se:.2f})\")\n",
    "\n",
    "# Visualization\n",
    "fig, axes = plt.subplots(2, 2, figsize=(12, 10))\n",
    "\n",
    "for idx, n in enumerate(sample_sizes):\n",
    "row = idx // 2\n",
    "col = idx % 2\n",
    "\n",
    "axes[row, col].hist(results[n], bins=50, alpha=0.7, color='blue', edgecolor='black')\n",
    "axes[row, col].set_title(f'Sample Size = {n}')\n",
    "axes[row, col].set_xlabel('Sample Mean')\n",
    "axes[row, col].set_ylabel('Frequency')\n",
    "axes[row, col].axvline(population.mean(), color='red', linestyle='--', label='Population Mean')\n",
    "axes[row, col].legend()\n",
    "axes[row, col].grid(True, alpha=0.3)\n",
    "\n",
    "plt.suptitle('Sampling Distribution of Means (Different Sample Sizes)', fontsize=14, fontweight='bold')\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "# Standard Error for different statistics\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"STANDARD ERROR FORMULAS\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "For sample mean: SE = σ / √n\n",
    "For sample proportion: SE = √(p(1-p) / n)\n",
    "For sample variance: More complex formula\n",
    "\"\"\")\n",
    "\n",
    "# Example: Proportion\n",
    "p = 0.4 # True population proportion\n",
    "n = 100\n",
    "se_proportion = np.sqrt(p * (1-p) / n)\n",
    "print(f\"\\nFor proportion: p={p}, n={n}\")\n",
    "print(f\"SE = √({p}×{1-p}/{n}) = {se_proportion:.4f}\")\n",
    "\n",
    "# Simulate\n",
    "sample_proportions = []\n",
    "for i in range(10000):\n",
    "sample = np.random.binomial(n, p) / n\n",
    "sample_proportions.append(sample)\n",
    "\n",
    "print(f\"Simulated SE: {np.std(sample_proportions):.4f}\")\n",
    "print(f\"Theoretical SE: {se_proportion:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83b6720c",
   "metadata": {},
   "source": [
    "### 2.2 The Central Limit Theorem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a22bbd2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print(\"=\"*50)\n",
    "print(\"CENTRAL LIMIT THEOREM (CLT)\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "THE CENTRAL LIMIT THEOREM STATES:\n",
    "If you take repeated samples from ANY distribution,\n",
    "the DISTRIBUTION OF THE SAMPLE MEANS approaches\n",
    "a NORMAL DISTRIBUTION as sample size increases.\n",
    "\n",
    "This is true regardless of the underlying distribution!\n",
    "\n",
    "Implications:\n",
    "- We can make inferences about populations using sample means\n",
    "- Normal distribution can approximate many distributions\n",
    "- Sample size matters: larger samples → more normal\n",
    "\"\"\")\n",
    "\n",
    "# Demonstrate CLT with different distributions\n",
    "fig, axes = plt.subplots(3, 3, figsize=(15, 12))\n",
    "\n",
    "# 1. UNIFORM DISTRIBUTION\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"1. UNIFORM DISTRIBUTION\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "# Generate population (uniform)\n",
    "population_uniform = np.random.uniform(0, 10, 50000)\n",
    "axes[0, 0].hist(population_uniform, bins=50, alpha=0.7, color='blue', edgecolor='black')\n",
    "axes[0, 0].set_title('Population (Uniform)')\n",
    "axes[0, 0].set_xlabel('Value')\n",
    "axes[0, 0].set_ylabel('Frequency')\n",
    "axes[0, 0].grid(True, alpha=0.3)\n",
    "\n",
    "# Sample means with n=10\n",
    "sample_means_10 = []\n",
    "for i in range(5000):\n",
    "sample = np.random.choice(population_uniform, 10)\n",
    "sample_means_10.append(sample.mean())\n",
    "\n",
    "axes[0, 1].hist(sample_means_10, bins=50, alpha=0.7, color='green', edgecolor='black')\n",
    "axes[0, 1].set_title('Sample Means (n=10)')\n",
    "axes[0, 1].set_xlabel('Sample Mean')\n",
    "axes[0, 1].set_ylabel('Frequency')\n",
    "axes[0, 1].grid(True, alpha=0.3)\n",
    "\n",
    "# Sample means with n=100\n",
    "sample_means_100 = []\n",
    "for i in range(5000):\n",
    "sample = np.random.choice(population_uniform, 100)\n",
    "sample_means_100.append(sample.mean())\n",
    "\n",
    "axes[0, 2].hist(sample_means_100, bins=50, alpha=0.7, color='red', edgecolor='black')\n",
    "axes[0, 2].set_title('Sample Means (n=100)')\n",
    "axes[0, 2].set_xlabel('Sample Mean')\n",
    "axes[0, 2].set_ylabel('Frequency')\n",
    "axes[0, 2].grid(True, alpha=0.3)\n",
    "\n",
    "# 2. EXPONENTIAL DISTRIBUTION (Right-skewed)\n",
    "print(\"\\n2. EXPONENTIAL DISTRIBUTION\")\n",
    "\n",
    "population_exp = np.random.exponential(2, 50000)\n",
    "axes[1, 0].hist(population_exp, bins=50, alpha=0.7, color='blue', edgecolor='black')\n",
    "axes[1, 0].set_title('Population (Exponential)')\n",
    "axes[1, 0].set_xlabel('Value')\n",
    "axes[1, 0].set_ylabel('Frequency')\n",
    "axes[1, 0].grid(True, alpha=0.3)\n",
    "\n",
    "sample_means_exp_10 = []\n",
    "for i in range(5000):\n",
    "sample = np.random.choice(population_exp, 10)\n",
    "sample_means_exp_10.append(sample.mean())\n",
    "\n",
    "axes[1, 1].hist(sample_means_exp_10, bins=50, alpha=0.7, color='green', edgecolor='black')\n",
    "axes[1, 1].set_title('Sample Means (n=10)')\n",
    "axes[1, 1].set_xlabel('Sample Mean')\n",
    "axes[1, 1].set_ylabel('Frequency')\n",
    "axes[1, 1].grid(True, alpha=0.3)\n",
    "\n",
    "sample_means_exp_100 = []\n",
    "for i in range(5000):\n",
    "sample = np.random.choice(population_exp, 100)\n",
    "sample_means_exp_100.append(sample.mean())\n",
    "\n",
    "axes[1, 2].hist(sample_means_exp_100, bins=50, alpha=0.7, color='red', edgecolor='black')\n",
    "axes[1, 2].set_title('Sample Means (n=100)')\n",
    "axes[1, 2].set_xlabel('Sample Mean')\n",
    "axes[1, 2].set_ylabel('Frequency')\n",
    "axes[1, 2].grid(True, alpha=0.3)\n",
    "\n",
    "# 3. BIMODAL DISTRIBUTION\n",
    "print(\"\\n3. BIMODAL DISTRIBUTION\")\n",
    "\n",
    "population_bimodal = np.concatenate([\n",
    "np.random.normal(20, 2, 25000),\n",
    "np.random.normal(80, 2, 25000)\n",
    "])\n",
    "\n",
    "axes[2, 0].hist(population_bimodal, bins=50, alpha=0.7, color='blue', edgecolor='black')\n",
    "axes[2, 0].set_title('Population (Bimodal)')\n",
    "axes[2, 0].set_xlabel('Value')\n",
    "axes[2, 0].set_ylabel('Frequency')\n",
    "axes[2, 0].grid(True, alpha=0.3)\n",
    "\n",
    "sample_means_bi_10 = []\n",
    "for i in range(5000):\n",
    "sample = np.random.choice(population_bimodal, 10)\n",
    "sample_means_bi_10.append(sample.mean())\n",
    "\n",
    "axes[2, 1].hist(sample_means_bi_10, bins=50, alpha=0.7, color='green', edgecolor='black')\n",
    "axes[2, 1].set_title('Sample Means (n=10)')\n",
    "axes[2, 1].set_xlabel('Sample Mean')\n",
    "axes[2, 1].set_ylabel('Frequency')\n",
    "axes[2, 1].grid(True, alpha=0.3)\n",
    "\n",
    "sample_means_bi_100 = []\n",
    "for i in range(5000):\n",
    "sample = np.random.choice(population_bimodal, 100)\n",
    "sample_means_bi_100.append(sample.mean())\n",
    "\n",
    "axes[2, 2].hist(sample_means_bi_100, bins=50, alpha=0.7, color='red', edgecolor='black')\n",
    "axes[2, 2].set_title('Sample Means (n=100)')\n",
    "axes[2, 2].set_xlabel('Sample Mean')\n",
    "axes[2, 2].set_ylabel('Frequency')\n",
    "axes[2, 2].grid(True, alpha=0.3)\n",
    "\n",
    "plt.suptitle('Central Limit Theorem: Sample Means Approach Normal Distribution',\n",
    "fontsize=14, fontweight='bold')\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"OBSERVATIONS\")\n",
    "print(\"=\"*50)\n",
    "print(\"As sample size increases, sample means distribution:\")\n",
    "print(\"1. Becomes more bell-shaped (normal)\")\n",
    "print(\"2. Becomes more concentrated (smaller variance)\")\n",
    "print(\"3. Approaches normal regardless of population shape\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78188722",
   "metadata": {},
   "source": [
    "## SESSION 3: Confidence Intervals, Bootstrapping, and Bayesian Thinking"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b4c20ed",
   "metadata": {},
   "source": [
    "### Duration: 2 hours"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3774bd05",
   "metadata": {},
   "source": [
    "### 3.1 Confidence Intervals and Point Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13e88390",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print(\"=\"*50)\n",
    "print(\"CONFIDENCE INTERVALS\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "WHAT IS A CONFIDENCE INTERVAL?\n",
    "- A range of values estimated to contain the true population parameter\n",
    "- Constructed from sample data\n",
    "- Associated with a confidence level (e.g., 95%)\n",
    "\n",
    "INTERPRETATION:\n",
    "If we repeated our sampling 100 times and created CI each time,\n",
    "95% of those intervals would contain the true population parameter.\n",
    "\n",
    "Common confidence levels: 90%, 95%, 99%\n",
    "\"\"\")\n",
    "\n",
    "# Example 1: Confidence interval for mean\n",
    "print(\"\\nEXAMPLE 1: CI FOR POPULATION MEAN\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "# Sample data\n",
    "sample_data = np.array([85, 92, 78, 88, 95, 82, 90, 87, 93, 80])\n",
    "print(f\"Sample: {sample_data}\")\n",
    "print(f\"n = {len(sample_data)}\")\n",
    "\n",
    "# Calculate statistics\n",
    "sample_mean = sample_data.mean()\n",
    "sample_std = sample_data.std(ddof=1) # Sample std dev\n",
    "se = sample_std / np.sqrt(len(sample_data))\n",
    "\n",
    "print(f\"Sample mean: {sample_mean:.2f}\")\n",
    "print(f\"Sample std dev: {sample_std:.2f}\")\n",
    "print(f\"Standard error: {se:.2f}\")\n",
    "\n",
    "# 95% Confidence interval using t-distribution\n",
    "confidence_level = 0.95\n",
    "alpha = 1 - confidence_level\n",
    "t_critical = stats.t.ppf(1 - alpha/2, df=len(sample_data)-1)\n",
    "margin_error = t_critical * se\n",
    "\n",
    "ci_lower = sample_mean - margin_error\n",
    "ci_upper = sample_mean + margin_error\n",
    "\n",
    "print(f\"\\nt-critical (df={len(sample_data)-1}): {t_critical:.3f}\")\n",
    "print(f\"Margin of error: {margin_error:.2f}\")\n",
    "print(f\"95% CI: [{ci_lower:.2f}, {ci_upper:.2f}]\")\n",
    "\n",
    "# Using scipy\n",
    "ci = stats.t.interval(0.95, len(sample_data)-1, loc=sample_mean, scale=se)\n",
    "print(f\"Using scipy: [{ci[0]:.2f}, {ci[1]:.2f}]\")\n",
    "\n",
    "# Example 2: Confidence interval for proportion\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"EXAMPLE 2: CI FOR PROPORTION\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "# Survey: 200 people, 120 prefer product A\n",
    "n = 200\n",
    "x = 120 # Number preferring A\n",
    "p_hat = x / n\n",
    "\n",
    "print(f\"Sample size: {n}\")\n",
    "print(f\"Number preferring: {x}\")\n",
    "print(f\"Sample proportion: {p_hat:.3f}\")\n",
    "\n",
    "# 95% CI for proportion\n",
    "se_prop = np.sqrt(p_hat * (1 - p_hat) / n)\n",
    "z_critical = stats.norm.ppf(0.975) # For 95%\n",
    "margin_error_prop = z_critical * se_prop\n",
    "\n",
    "ci_lower_prop = p_hat - margin_error_prop\n",
    "ci_upper_prop = p_hat + margin_error_prop\n",
    "\n",
    "print(f\"Standard error: {se_prop:.4f}\")\n",
    "print(f\"z-critical: {z_critical:.3f}\")\n",
    "print(f\"Margin of error: {margin_error_prop:.3f}\")\n",
    "print(f\"95% CI: [{ci_lower_prop:.3f}, {ci_upper_prop:.3f}]\")\n",
    "print(f\"95% CI (percentage): [{ci_lower_prop*100:.1f}%, {ci_upper_prop*100:.1f}%]\")\n",
    "\n",
    "# Effect of confidence level\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"EFFECT OF CONFIDENCE LEVEL\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "confidence_levels = [0.90, 0.95, 0.99]\n",
    "for conf_level in confidence_levels:\n",
    "alpha = 1 - conf_level\n",
    "t_crit = stats.t.ppf(1 - alpha/2, df=len(sample_data)-1)\n",
    "me = t_crit * se\n",
    "ci_l = sample_mean - me\n",
    "ci_u = sample_mean + me\n",
    "print(f\"{conf_level*100:.0f}% CI: [{ci_l:.2f}, {ci_u:.2f}], Width: {ci_u-ci_l:.2f}\")\n",
    "\n",
    "print(\"\\nObservation: Higher confidence level → Wider interval\")\n",
    "\n",
    "# Visualization\n",
    "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
    "\n",
    "# Plot 1: Sample means and CI\n",
    "np.random.seed(42)\n",
    "sample_means = []\n",
    "cis_lower = []\n",
    "cis_upper = []\n",
    "\n",
    "true_mean = 100\n",
    "for i in range(50):\n",
    "sample = np.random.normal(true_mean, 15, 30)\n",
    "m = sample.mean()\n",
    "s = sample.std(ddof=1)\n",
    "se_i = s / np.sqrt(30)\n",
    "t_crit = stats.t.ppf(0.975, df=29)\n",
    "\n",
    "sample_means.append(m)\n",
    "cis_lower.append(m - t_crit * se_i)\n",
    "cis_upper.append(m + t_crit * se_i)\n",
    "\n",
    "colors = ['red' if true_mean < cis_lower[i] or true_mean > cis_upper[i] else 'blue'\n",
    "for i in range(50)]\n",
    "\n",
    "for i in range(50):\n",
    "axes[0].plot([cis_lower[i], cis_upper[i]], [i, i], color=colors[i], linewidth=2)\n",
    "axes[0].scatter(sample_means[i], i, color=colors[i], s=50, zorder=5)\n",
    "\n",
    "axes[0].axvline(true_mean, color='green', linestyle='--', linewidth=2, label='True Mean')\n",
    "axes[0].set_xlabel('Value')\n",
    "axes[0].set_ylabel('Sample Number')\n",
    "axes[0].set_title('50 Confidence Intervals (95%)\n",
    "(Blue include true mean, Red do not)')\n",
    "axes[0].legend()\n",
    "axes[0].grid(True, alpha=0.3)\n",
    "\n",
    "# Plot 2: Effect of sample size\n",
    "sample_sizes = [10, 30, 50, 100, 200]\n",
    "margin_errors = []\n",
    "\n",
    "for n in sample_sizes:\n",
    "se_n = 15 / np.sqrt(n)\n",
    "t_crit = stats.t.ppf(0.975, df=n-1)\n",
    "me = t_crit * se_n\n",
    "margin_errors.append(me)\n",
    "\n",
    "axes[1].plot(sample_sizes, margin_errors, marker='o', linewidth=2, markersize=8, color='blue')\n",
    "axes[1].set_xlabel('Sample Size')\n",
    "axes[1].set_ylabel('Margin of Error')\n",
    "axes[1].set_title('Effect of Sample Size on Margin of Error\n",
    "(95% CI)')\n",
    "axes[1].grid(True, alpha=0.3)\n",
    "axes[1].set_xscale('log')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53e30178",
   "metadata": {},
   "source": [
    "### 3.2 Bootstrapping for Statistical Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3bbff63",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print(\"=\"*50)\n",
    "print(\"BOOTSTRAPPING\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "WHAT IS BOOTSTRAPPING?\n",
    "- Resampling method: repeatedly sample from data with replacement\n",
    "- Estimates distribution of a statistic without assumptions\n",
    "- Very powerful for complex statistics where theoretical distributions unknown\n",
    "\n",
    "STEPS:\n",
    "1. Draw random sample WITH REPLACEMENT from original sample\n",
    "2. Calculate statistic of interest on resampled data\n",
    "3. Repeat steps 1-2 many times (1000-10000)\n",
    "4. Use distribution of results to estimate CI and SE\n",
    "\"\"\")\n",
    "\n",
    "# Original data\n",
    "original_data = np.array([85, 92, 78, 88, 95, 82, 90, 87, 93, 80, 86, 91])\n",
    "print(f\"\\nOriginal data: {original_data}\")\n",
    "print(f\"Original mean: {original_data.mean():.2f}\")\n",
    "print(f\"Original median: {np.median(original_data):.2f}\")\n",
    "\n",
    "# Bootstrap for mean\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"BOOTSTRAP: MEAN\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "n_bootstrap = 10000\n",
    "bootstrap_means = []\n",
    "\n",
    "np.random.seed(42)\n",
    "for i in range(n_bootstrap):\n",
    "# Resample with replacement\n",
    "bootstrap_sample = np.random.choice(original_data, len(original_data), replace=True)\n",
    "bootstrap_means.append(bootstrap_sample.mean())\n",
    "\n",
    "bootstrap_means = np.array(bootstrap_means)\n",
    "\n",
    "print(f\"\\nBootstrap distribution statistics:\")\n",
    "print(f\"Mean of bootstrap means: {bootstrap_means.mean():.2f}\")\n",
    "print(f\"Std of bootstrap means (SE): {bootstrap_means.std():.2f}\")\n",
    "\n",
    "# 95% CI using percentile method\n",
    "ci_lower = np.percentile(bootstrap_means, 2.5)\n",
    "ci_upper = np.percentile(bootstrap_means, 97.5)\n",
    "print(f\"95% CI (percentile method): [{ci_lower:.2f}, {ci_upper:.2f}]\")\n",
    "\n",
    "# Bootstrap for median\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"BOOTSTRAP: MEDIAN\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "bootstrap_medians = []\n",
    "for i in range(n_bootstrap):\n",
    "bootstrap_sample = np.random.choice(original_data, len(original_data), replace=True)\n",
    "bootstrap_medians.append(np.median(bootstrap_sample))\n",
    "\n",
    "bootstrap_medians = np.array(bootstrap_medians)\n",
    "\n",
    "print(f\"\\nBootstrap median statistics:\")\n",
    "print(f\"Mean of bootstrap medians: {bootstrap_medians.mean():.2f}\")\n",
    "print(f\"Std of bootstrap medians (SE): {bootstrap_medians.std():.2f}\")\n",
    "\n",
    "ci_lower_median = np.percentile(bootstrap_medians, 2.5)\n",
    "ci_upper_median = np.percentile(bootstrap_medians, 97.5)\n",
    "print(f\"95% CI for median: [{ci_lower_median:.2f}, {ci_upper_median:.2f}]\")\n",
    "\n",
    "# Bootstrap for correlation\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"BOOTSTRAP: CORRELATION\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "# Create paired data\n",
    "np.random.seed(42)\n",
    "x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])\n",
    "y = 2 * x + np.random.normal(0, 2, 10)\n",
    "\n",
    "original_corr = np.corrcoef(x, y)[0, 1]\n",
    "print(f\"Original correlation: {original_corr:.3f}\")\n",
    "\n",
    "bootstrap_corrs = []\n",
    "for i in range(n_bootstrap):\n",
    "indices = np.random.choice(len(x), len(x), replace=True)\n",
    "x_boot = x[indices]\n",
    "y_boot = y[indices]\n",
    "corr_boot = np.corrcoef(x_boot, y_boot)[0, 1]\n",
    "bootstrap_corrs.append(corr_boot)\n",
    "\n",
    "bootstrap_corrs = np.array(bootstrap_corrs)\n",
    "\n",
    "print(f\"Bootstrap correlation mean: {bootstrap_corrs.mean():.3f}\")\n",
    "print(f\"Bootstrap correlation SE: {bootstrap_corrs.std():.3f}\")\n",
    "\n",
    "ci_lower_corr = np.percentile(bootstrap_corrs, 2.5)\n",
    "ci_upper_corr = np.percentile(bootstrap_corrs, 97.5)\n",
    "print(f\"95% CI for correlation: [{ci_lower_corr:.3f}, {ci_upper_corr:.3f}]\")\n",
    "\n",
    "# Visualization\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 5))\n",
    "\n",
    "# Bootstrap means\n",
    "axes[0].hist(bootstrap_means, bins=50, alpha=0.7, color='blue', edgecolor='black')\n",
    "axes[0].axvline(original_data.mean(), color='red', linestyle='--', linewidth=2, label='Original Mean')\n",
    "axes[0].axvline(ci_lower, color='green', linestyle='--', linewidth=2, label='95% CI')\n",
    "axes[0].axvline(ci_upper, color='green', linestyle='--', linewidth=2)\n",
    "axes[0].set_title('Bootstrap Distribution of Mean')\n",
    "axes[0].set_xlabel('Sample Mean')\n",
    "axes[0].set_ylabel('Frequency')\n",
    "axes[0].legend()\n",
    "axes[0].grid(True, alpha=0.3)\n",
    "\n",
    "# Bootstrap medians\n",
    "axes[1].hist(bootstrap_medians, bins=50, alpha=0.7, color='green', edgecolor='black')\n",
    "axes[1].axvline(np.median(original_data), color='red', linestyle='--', linewidth=2, label='Original Median')\n",
    "axes[1].axvline(ci_lower_median, color='orange', linestyle='--', linewidth=2, label='95% CI')\n",
    "axes[1].axvline(ci_upper_median, color='orange', linestyle='--', linewidth=2)\n",
    "axes[1].set_title('Bootstrap Distribution of Median')\n",
    "axes[1].set_xlabel('Sample Median')\n",
    "axes[1].set_ylabel('Frequency')\n",
    "axes[1].legend()\n",
    "axes[1].grid(True, alpha=0.3)\n",
    "\n",
    "# Bootstrap correlations\n",
    "axes[2].hist(bootstrap_corrs, bins=50, alpha=0.7, color='purple', edgecolor='black')\n",
    "axes[2].axvline(original_corr, color='red', linestyle='--', linewidth=2, label='Original Correlation')\n",
    "axes[2].axvline(ci_lower_corr, color='orange', linestyle='--', linewidth=2, label='95% CI')\n",
    "axes[2].axvline(ci_upper_corr, color='orange', linestyle='--', linewidth=2)\n",
    "axes[2].set_title('Bootstrap Distribution of Correlation')\n",
    "axes[2].set_xlabel('Sample Correlation')\n",
    "axes[2].set_ylabel('Frequency')\n",
    "axes[2].legend()\n",
    "axes[2].grid(True, alpha=0.3)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "863554fd",
   "metadata": {},
   "source": [
    "### 3.3 Bayesian Thinking and Prior Knowledge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "759d3941",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print(\"=\"*50)\n",
    "print(\"BAYESIAN THINKING\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "BAYES' THEOREM:\n",
    "P(A|B) = P(B|A) × P(A) / P(B)\n",
    "\n",
    "COMPONENTS:\n",
    "- P(A|B) = Posterior (what we want to know)\n",
    "- P(B|A) = Likelihood (probability of evidence given hypothesis)\n",
    "- P(A) = Prior (what we believed before data)\n",
    "- P(B) = Marginal likelihood (normalizing factor)\n",
    "\n",
    "KEY INSIGHT:\n",
    "Bayesian approach updates beliefs as we get new data\n",
    "\"\"\")\n",
    "\n",
    "# Example 1: Disease testing revisited\n",
    "print(\"\\nEXAMPLE 1: MEDICAL TESTING\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "# Prior: 1% of population has disease\n",
    "prior_disease = 0.01\n",
    "prior_no_disease = 1 - prior_disease\n",
    "\n",
    "# Test accuracy\n",
    "sensitivity = 0.95 # P(+|disease)\n",
    "specificity = 0.90 # P(-|no disease)\n",
    "\n",
    "# What we observe: positive test\n",
    "# P(+|disease) × P(disease) = likelihood × prior\n",
    "prob_positive_given_disease = sensitivity * prior_disease\n",
    "prob_positive_given_no_disease = (1 - specificity) * prior_no_disease\n",
    "\n",
    "# P(+) = total probability of positive\n",
    "prob_positive_total = prob_positive_given_disease + prob_positive_given_no_disease\n",
    "\n",
    "# Posterior: P(disease|+)\n",
    "posterior_disease = prob_positive_given_disease / prob_positive_total\n",
    "\n",
    "print(f\"Prior P(disease) = {prior_disease:.3f}\")\n",
    "print(f\"Sensitivity P(+|disease) = {sensitivity:.3f}\")\n",
    "print(f\"Specificity P(-|no disease) = {specificity:.3f}\")\n",
    "print(f\"False positive rate = {1-specificity:.3f}\")\n",
    "print(f\"\\nAfter positive test:\")\n",
    "print(f\"Posterior P(disease|+) = {posterior_disease:.3f} or {posterior_disease*100:.1f}%\")\n",
    "print(f\"\\nInterpretation: Even with positive test, only {posterior_disease*100:.1f}% chance of disease!\")\n",
    "\n",
    "# Example 2: Coin fairness\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"EXAMPLE 2: IS COIN FAIR?\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "# Prior: assume coin is fair (50-50)\n",
    "# Data: flip 10 times, get 8 heads\n",
    "\n",
    "print(f\"Prior belief: coin is fair (P(heads)=0.5)\")\n",
    "print(f\"Data: 10 flips, 8 heads\")\n",
    "\n",
    "# Calculate likelihood for different values of p\n",
    "p_values = np.linspace(0, 1, 100)\n",
    "data = 8 # heads\n",
    "trials = 10\n",
    "\n",
    "# Likelihood (binomial)\n",
    "likelihood = stats.binom.pmf(data, trials, p_values)\n",
    "\n",
    "# Uniform prior\n",
    "prior = np.ones_like(p_values)\n",
    "\n",
    "# Posterior (proportional to likelihood × prior)\n",
    "posterior = likelihood * prior\n",
    "posterior = posterior / posterior.sum() # Normalize\n",
    "\n",
    "# Plot\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.plot(p_values, likelihood, label='Likelihood (Data: 8H in 10 flips)', linewidth=2)\n",
    "plt.plot(p_values, prior, label='Prior (uniform)', linewidth=2)\n",
    "plt.plot(p_values, posterior, label='Posterior', linewidth=3)\n",
    "\n",
    "plt.axvline(0.5, color='red', linestyle='--', label='Fair coin (p=0.5)')\n",
    "plt.xlabel('Probability of Heads (p)')\n",
    "plt.ylabel('Probability Density')\n",
    "plt.title('Bayesian Inference: Is the Coin Fair?')\n",
    "plt.legend(fontsize=11)\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "# Calculate posterior mean\n",
    "posterior_mean = np.average(p_values, weights=posterior)\n",
    "print(f\"\\nPosterior mean (best estimate): {posterior_mean:.3f}\")\n",
    "print(f\"Our estimate shifted from 0.5 (prior) to {posterior_mean:.3f} (posterior)\")\n",
    "\n",
    "# Example 3: Updating beliefs\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"EXAMPLE 3: SEQUENTIAL LEARNING\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(f\"Starting belief: coin is fair (p=0.5)\")\n",
    "print(f\"\\nSequence of events:\")\n",
    "\n",
    "p = 0.8 # True probability\n",
    "prior_mean = 0.5\n",
    "\n",
    "observations = [1, 1, 0, 1, 1, 1, 0, 1, 1, 1] # 1=heads, 0=tails\n",
    "\n",
    "for i, obs in enumerate(observations, 1):\n",
    "# Count heads so far\n",
    "heads_count = sum(observations[:i])\n",
    "\n",
    "# Update using Beta-Binomial model\n",
    "# Beta distribution for prior on p\n",
    "alpha = 1 + heads_count\n",
    "beta = 1 + (i - heads_count)\n",
    "\n",
    "# Posterior mean\n",
    "posterior_mean_beta = alpha / (alpha + beta)\n",
    "\n",
    "outcome = \"H\" if obs == 1 else \"T\"\n",
    "print(f\"Flip {i:2d}: {outcome} → Heads: {heads_count}, Posterior mean: {posterior_mean_beta:.3f}\")\n",
    "\n",
    "print(f\"\\nFinal posterior mean: {posterior_mean_beta:.3f}\")\n",
    "print(f\"(True value was p={p})\")\n",
    "\n",
    "# Compare Frequentist vs Bayesian\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"FREQUENTIST vs BAYESIAN\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(f\"\"\"\n",
    "FREQUENTIST:\n",
    "- Probability is long-run frequency\n",
    "- Parameter is fixed, unknown\n",
    "- Uses only the data\n",
    "- P-values and confidence intervals\n",
    "\n",
    "BAYESIAN:\n",
    "- Probability is subjective degree of belief\n",
    "- Parameter has a distribution\n",
    "- Combines data with prior knowledge\n",
    "- Posterior distributions and credible intervals\n",
    "\n",
    "WHEN TO USE:\n",
    "- Frequentist: When you have lots of data, want standard tests\n",
    "- Bayesian: When you want to incorporate prior knowledge, small samples\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b11f5541",
   "metadata": {},
   "source": [
    "## Week 6 Summary\n",
    "By completing Week 6, you have learned:\n",
    "- Probability fundamentals: sample space, events, basic rules\n",
    "- Addition rule: P(A or B)\n",
    "- Multiplication rule: P(A and B)\n",
    "- Conditional probability and Bayes&#x27; theorem\n",
    "- Discrete distributions: Bernoulli, Binomial, Poisson\n",
    "- Continuous distributions: Normal, Exponential, Uniform\n",
    "- Probability mass function (PMF) and probability density function (PDF)\n",
    "- Working with scipy.stats for distributions\n",
    "- Sampling distributions and standard error\n",
    "- Central Limit Theorem: why sampling distributions approach normal\n",
    "- Confidence intervals for means and proportions\n",
    "- Margin of error and confidence levels\n",
    "- Bootstrapping: resampling with replacement\n",
    "- Bootstrapping for complex statistics\n",
    "- Bayesian thinking: priors, likelihoods, posteriors\n",
    "- Comparing Frequentist vs Bayesian approaches"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6dd01317",
   "metadata": {},
   "source": [
    "## Week 6 Assignments"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd90aef0",
   "metadata": {},
   "source": [
    "### Assignment 1: Probability and Distributions\n",
    "Work with probability distributions:\n",
    "- Calculate probabilities for various discrete distributions\n",
    "- Solve conditional probability and Bayes theorem problems\n",
    "- Plot and compare probability distributions\n",
    "- Demonstrate Central Limit Theorem with simulations\n",
    "- Calculate probabilities for real-world scenarios\n",
    "- Create visualizations showing different distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11696707",
   "metadata": {},
   "source": [
    "### Assignment 2: Confidence Intervals and Bootstrapping\n",
    "Perform interval estimation:\n",
    "- Calculate confidence intervals for different statistics\n",
    "- Interpret confidence levels correctly\n",
    "- Perform bootstrap resampling to estimate CIs\n",
    "- Compare bootstrap CIs with theoretical CIs\n",
    "- Explore effect of sample size on margin of error\n",
    "- Create visualizations of confidence intervals"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcecf18f",
   "metadata": {},
   "source": [
    "### Assignment 3: Comprehensive Statistical Analysis\n",
    "Apply all concepts to real data:\n",
    "- Load a dataset and perform comprehensive statistical analysis\n",
    "- Calculate point estimates and confidence intervals\n",
    "- Use bootstrap for complex statistics\n",
    "- Demonstrate understanding of Central Limit Theorem\n",
    "- Interpret results in business context\n",
    "- Create professional report with visualizations\n",
    "- Compare results using both frequentist and Bayesian approaches"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa03e7fb",
   "metadata": {},
   "source": [
    "## Recommended Practice Problems\n",
    "- Solve probability problems using Bayes theorem\n",
    "- Calculate probabilities for various distributions\n",
    "- Demonstrate CLT with different sample sizes\n",
    "- Compare theoretical and simulated distributions\n",
    "- Interpret confidence intervals correctly\n",
    "- Perform bootstrapping on different statistics\n",
    "- Solve problems combining multiple probability distributions\n",
    "- Create simulations to understand statistical concepts\n",
    "- Compare CI methods (theoretical vs bootstrap)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6d4a79e",
   "metadata": {},
   "source": [
    "## Additional Resources"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0afd8c3",
   "metadata": {},
   "source": [
    "### Books\n",
    "- Chapter 6: Probability and Statistics - &quot;Python for Data Analysis&quot; by Wes McKinney\n",
    "- Statistical Rethinking by Richard McElreath - Bayesian approach\n",
    "- Think Stats by Allen Downey - Free online"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05e4c230",
   "metadata": {},
   "source": [
    "### Online Documentation\n",
    "- scipy.stats documentation: https://docs.scipy.org/doc/scipy/reference/stats.html\n",
    "- NumPy random: https://numpy.org/doc/stable/reference/random/index.html\n",
    "- Bayesian Inference tutorial: https://en.wikipedia.org/wiki/Bayesian_inference\n",
    "- Central Limit Theorem visualizations"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}