{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "colab": {
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "cells": [
  {
   "id": "m02-l1-0",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lesson 2.1 · Data Generating Processes and Distributions\n",
    "\n",
    "*Accelerated Foundation in Applied Causal Inference*\n",
    "\n",
    "**By the end of this lesson you can:**\n",
    "\n",
    "- Define a data generating process (DGP) and explain why simulation replaces proofs in this course\n",
    "- Simulate continuous and discrete variables with NumPy and scipy.stats\n",
    "- Distinguish a population parameter from a sample statistic\n",
    "- Summarise and visualise a simulated cohort\n",
    "\n",
    "> Read the full explanations on the course website, then run each cell here in order (**Runtime → Run all** works on a fresh runtime)."
   ]
  },
  {
   "id": "m02-l1-1",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why we simulate: playing god with data"
   ]
  },
  {
   "id": "m02-l1-2",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup and the random number generator"
   ]
  },
  {
   "id": "m02-l1-3",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "!pip install -q numpy scipy polars matplotlib seaborn\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import polars as pl\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "rng = np.random.default_rng(101)     # a seeded generator: same numbers every run\n",
    "sns.set_theme(style=\"whitegrid\")\n",
    "print(rng.normal(size=3))"
   ]
  },
  {
   "id": "m02-l1-4",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building a public-health DGP"
   ]
  },
  {
   "id": "m02-l1-5",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question it answers:** What does a population of 50,000 patients look like if age is normal and 25% have a comorbidity?"
   ]
  },
  {
   "id": "m02-l1-6",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "POPULATION_SIZE = 50_000\n",
    "TRUE_MEAN_AGE = 45.0\n",
    "TRUE_SD_AGE = 12.0\n",
    "TRUE_COMORBIDITY_RATE = 0.25\n",
    "\n",
    "age = np.clip(rng.normal(TRUE_MEAN_AGE, TRUE_SD_AGE, POPULATION_SIZE), 18, 90)\n",
    "comorbidity = rng.binomial(n=1, p=TRUE_COMORBIDITY_RATE, size=POPULATION_SIZE)\n",
    "visits = rng.poisson(lam=2 + 2 * comorbidity)          # GP visits per year\n",
    "\n",
    "population = pl.DataFrame({\n",
    "    \"patient_id\": np.arange(POPULATION_SIZE),\n",
    "    \"age\": age,\n",
    "    \"has_comorbidity\": comorbidity,\n",
    "    \"gp_visits\": visits,\n",
    "})\n",
    "print(population.head())\n",
    "print(population.select(\n",
    "    pl.col(\"age\").mean().round(2).alias(\"mean_age\"),\n",
    "    pl.col(\"has_comorbidity\").mean().round(3).alias(\"comorbidity_rate\"),\n",
    "    pl.col(\"gp_visits\").mean().round(2).alias(\"mean_visits\"),\n",
    "))"
   ]
  },
  {
   "id": "m02-l1-7",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(14, 4))\n",
    "sns.histplot(population[\"age\"].to_numpy(), bins=40, ax=axes[0]).set(title=\"Age (normal, clipped)\")\n",
    "sns.countplot(x=population[\"has_comorbidity\"].to_numpy(), ax=axes[1]).set(title=\"Comorbidity (Bernoulli)\")\n",
    "sns.countplot(x=population[\"gp_visits\"].to_numpy(), ax=axes[2]).set(title=\"GP visits (Poisson)\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "id": "m02-l1-8",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distributions you will meet all course long"
   ]
  },
  {
   "id": "m02-l1-9",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question it answers:** What share of patients is older than 65, in theory and in the simulation?"
   ]
  },
  {
   "id": "m02-l1-10",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "theory = 1 - stats.norm.cdf(65, loc=TRUE_MEAN_AGE, scale=TRUE_SD_AGE)\n",
    "simulated = (population[\"age\"] > 65).mean()\n",
    "print(f\"P(age > 65): theory {theory:.4f} | simulation {simulated:.4f}\")"
   ]
  },
  {
   "id": "m02-l1-11",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: an income distribution"
   ]
  },
  {
   "id": "m02-l1-12",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ✍️ Exercise\n",
    "\n",
    "Simulate annual incomes for the 50,000 patients from a log-normal distribution with `mean=10.8` and\n",
    "`sigma=0.6` (on the log scale). Report the mean and the median income. Why are they so different? Which\n",
    "would you report in a policy memo about a \"typical\" household?"
   ]
  },
  {
   "id": "m02-l1-13",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Your answer here"
   ]
  },
  {
   "id": "m02-l1-14",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checklist and recap"
   ]
  }
 ]
}