{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "iid.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# independent and identically distributed (i.i.d.) \u2014 Python demo\n\nNumerical companion to the entry [independent and identically distributed (i.i.d.)](https://dictionaryofml.org/terms/iid.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nThe entry asks whether temperature measurements recorded at Krems in 2024 can be modeled as i.i.d. random variables, and answers it by testing the two halves of the definition separately on the record itself. Self-contained (numpy/matplotlib only), fixed seed; the measurements are fetched from the public GeoSphere Austria archive.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/iid.py`](https://dictionaryofml.org/terms/iid.py); CC BY 4.0."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# Notebook shim: the script resolves output paths relative to __file__,\n# which a notebook kernel does not define; everything lands in the\n# working directory instead.\nimport os\n__file__ = os.path.join(os.getcwd(), \"iid.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\niid.py \u2014 numerical companion to the glossary entry 'independent and\nidentically distributed (i.i.d.)'.\n\nThe entry asks whether temperature measurements recorded at Krems in 2024 can\nbe modeled as i.i.d. random variables, and answers it by testing the two halves\nof the definition separately on the record itself. Self-contained\n(numpy/matplotlib only), fixed seed; the measurements are fetched from the\npublic GeoSphere Austria archive.\n\nBlocks\n------\n[B-months]      Air temperature every ten minutes for January, May, August and\n                November, one dataset per month. The levels differ across the\n                months and the daily cycle repeats inside each of them, which\n                is what the rest of the entry measures.\n[B-data]        366 daily maximum temperatures at Krems an der Donau\n                (station 3805) for 2024.\n[B-identical]   The identically distributed half fails: the monthly averages\n                run from January to July, so the distribution of the\n                temperature depends on which day is read.\n[B-independent] The independent half fails too: the event \"above 25 degrees\"\n                has probability 0.265 per day, so under independence two\n                consecutive days would both exceed it with probability 0.070;\n                the record does it three times as often. The lag-one\n                correlation says the same.\n[B-shuffle]     A random permutation of the same 366 numbers leaves the\n                collection of values untouched and makes the product rule\n                hold, which is what independence is about. Writes the two\n                panels of the entry's second figure.\n[B-deseason]    Subtracting the seasonal average repairs the identically\n                distributed half only: the lag-one correlation of what is\n                left is still 0.67.\n[B-verify]      The two standard methods applied to periods of growing length\n                (January, January to June, the whole year): the largest gap\n                between the empirical CDFs of two 15-day blocks of the period\n                against the Kolmogorov-Smirnov threshold, and the lag-one\n                correlation calibrated by a permutation test. Writes the\n                panels of the entry's third figure.\n[B-optimal]     The same statistic on one month and on one year: the\n                correlation of consecutive values over random reorderings,\n                which tightens\n                as the collection grows. Writes the panels of the entry's\n                last figure.\n\nOutputs\n-------\niid_temps.csv            : date, daily maximum temperature, 366 days of 2024\niid_month_<mon>.csv      : day, temperature every ten minutes, four months\niid_lag_record.csv       : today, tomorrow -- the record, for the left panel\niid_lag_shuffled.csv     : today, tomorrow -- one permutation, right panel\niid_cdf_<period>_<block>.csv : temperature, fraction -- empirical CDFs\niid_perm_<period>.csv    : center, fraction -- the correlation of consecutive\n                           values over random reorderings\niid.png                  : preview (checking only)\n\nData generated by pythondemos/iid.py.\n\"\"\"\n\nimport json\nimport math\nimport urllib.request\nfrom pathlib import Path\n\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nOUT_DIR = Path(__file__).parent\nARCHIVE = \"https://dataset.api.hub.geosphere.at/v1/station/historical/\"\nSTATION = 3805\n\nreport = []\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")\n\n\ndef fetch(resource, parameter, start, end):\n    \"\"\"One parameter of one station from the GeoSphere Austria archive.\"\"\"\n    url = (f\"{ARCHIVE}{resource}?parameters={parameter}\"\n           f\"&station_ids={STATION}&start={start}&end={end}\")\n    with urllib.request.urlopen(url, timeout=180) as resp:\n        payload = json.load(resp)\n    values = payload[\"features\"][0][\"properties\"][\"parameters\"]\n    key = list(values.keys())[0]\n    return (np.array(values[key][\"data\"], dtype=float),\n            [t[:16] for t in payload[\"timestamps\"]])\n\n\ndef lag_one_correlation(x):\n    \"\"\"Pearson correlation coefficient between consecutive entries of x.\"\"\"\n    return float(np.corrcoef(x[:-1], x[1:])[0, 1])\n\n\ndef ks_statistic(a, b):\n    \"\"\"Largest gap between the empirical CDFs of two collections.\"\"\"\n    grid = np.sort(np.concatenate([a, b]))\n    fa = np.searchsorted(np.sort(a), grid, side=\"right\") / len(a)\n    fb = np.searchsorted(np.sort(b), grid, side=\"right\") / len(b)\n    return float(np.abs(fa - fb).max())\n\n\ndef ks_threshold(size, level):\n    \"\"\"Two-sample Kolmogorov-Smirnov critical value for two blocks of size.\"\"\"\n    return math.sqrt(-0.5 * math.log(level / 2)) * math.sqrt(2 / size)\n\n\ndef empirical_cdf(x):\n    \"\"\"Sorted values and the fraction of the collection at or below each.\"\"\"\n    values = np.sort(x)\n    return values, np.arange(1, len(values) + 1) / len(values)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-months]** Air temperature every ten minutes for January, May, August and November, one dataset per month. The levels differ across the months and the daily cycle repeats inside each of them, which is what the rest of the entry measures."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "MONTHS = [(1, \"jan\", \"January\"), (5, \"may\", \"May\"),\n          (8, \"aug\", \"August\"), (11, \"nov\", \"November\")]\nfine = {}\nfor number, tag, title in MONTHS:\n    last = 31 if number in (1, 8) else 30\n    values, stamps = fetch(\"klima-v2-10min\", \"TL\",\n                           f\"2024-{number:02d}-01T00:00\",\n                           f\"2024-{number:02d}-{last:02d}T23:50\")\n    day = np.array([int(s[8:10]) + int(s[11:13]) / 24 + int(s[14:16]) / 1440\n                    for s in stamps])\n    fine[tag] = (day, values, title)\n    with open(OUT_DIR / f\"iid_month_{tag}.csv\", \"w\") as f:\n        f.write(\"day,temp\\n\")\n        for a, b in zip(day, values):\n            f.write(f\"{a:.4f},{b:g}\\n\")\n    print(f\"[B-months] {title}: {len(values)} measurements, mean \"\n          f\"{values.mean():.2f} deg, lag-one correlation \"\n          f\"{lag_one_correlation(values):.3f}\")\ncheck(\"[B-months] every month is measured every ten minutes\",\n      all(len(fine[tag][1]) in (4320, 4464) for _, tag, _ in MONTHS))\ncheck(\"[B-months] the monthly levels differ by more than twenty degrees\",\n      fine[\"aug\"][1].mean() - fine[\"jan\"][1].mean() > 20.0)\ncheck(\"[B-months] consecutive measurements are nearly equal in every month\",\n      all(lag_one_correlation(fine[tag][1]) > 0.98 for _, tag, _ in MONTHS))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-data]** 366 daily maximum temperatures at Krems an der Donau (station 3805) for 2024."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "temp, stamps = fetch(\"klima-v2-1d\", \"tlmax\", \"2024-01-01\", \"2024-12-31\")\nmonth = np.array([int(s[5:7]) for s in stamps])\nwith open(OUT_DIR / \"iid_temps.csv\", \"w\") as f:\n    f.write(\"date,tlmax\\n\")\n    for day, value in zip(stamps, temp):\n        f.write(f\"{day[:10]},{value:g}\\n\")\nprint(f\"[B-data] {len(temp)} daily maxima, mean {temp.mean():.2f} deg, \"\n      f\"variance {temp.var(ddof=1):.2f} deg squared\")\ncheck(\"[B-data] the year has all 366 days\", len(temp) == 366)\ncheck(\"[B-data] the record matches the archive (Feb 1: 10.4 deg)\",\n      stamps[31][:10] == \"2024-02-01\" and np.isclose(temp[31], 10.4))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-identical]** The identically distributed half fails: the monthly averages run from January to July, so the distribution of the temperature depends on which day is read."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "jan, jul = temp[month == 1], temp[month == 7]\nprint(f\"[B-identical] January averages {jan.mean():.2f} deg and July \"\n      f\"{jul.mean():.2f} deg\")\ncheck(\"[B-identical] the July average exceeds the January average by more \"\n      \"than twenty degrees\", jul.mean() - jan.mean() > 20.0)\ncheck(\"[B-identical] the gap is large against the spread within a month\",\n      jul.mean() - jan.mean() > 4.0 * max(jan.std(ddof=1), jul.std(ddof=1)))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-independent]** The independent half fails too: the event \"above 25 degrees\" has probability 0.265 per day, so under independence two consecutive days would both exceed it with probability 0.070; the record does it three times as often. The lag-one correlation says the same."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "THRESHOLD = 25.0\nwarm = temp > THRESHOLD\np_single = warm.mean()\np_pair = (warm[:-1] & warm[1:]).mean()\nprint(f\"[B-independent] a day above {THRESHOLD:g} deg has frequency \"\n      f\"{p_single:.3f}; two in a row {p_pair:.3f} against the product \"\n      f\"{p_single ** 2:.3f}, a factor of {p_pair / p_single ** 2:.2f}\")\nprint(f\"[B-independent] lag-one correlation of the record \"\n      f\"{lag_one_correlation(temp):.3f}\")\ncheck(\"[B-independent] consecutive warm days are at least three times as \"\n      \"frequent as the product rule allows\", p_pair > 3.0 * p_single ** 2)\ncheck(\"[B-independent] the lag-one correlation is above 0.9\",\n      lag_one_correlation(temp) > 0.9)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-shuffle]** A random permutation of the same 366 numbers leaves the collection of values untouched and makes the product rule hold, which is what independence is about. Writes the two panels of the entry's second figure."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "rng = np.random.default_rng(0)\nshuffled = rng.permutation(temp)\nwarm_s = shuffled > THRESHOLD\np_pair_s = (warm_s[:-1] & warm_s[1:]).mean()\nprint(f\"[B-shuffle] after permuting: two warm days in a row {p_pair_s:.3f} \"\n      f\"against the product {p_single ** 2:.3f}; lag-one correlation \"\n      f\"{lag_one_correlation(shuffled):.3f}\")\ncheck(\"[B-shuffle] permuting leaves the collection of values unchanged\",\n      np.allclose(np.sort(shuffled), np.sort(temp)))\ncheck(\"[B-shuffle] the product rule now holds to within a tenth\",\n      abs(p_pair_s - p_single ** 2) < 0.1 * p_single ** 2)\ncheck(\"[B-shuffle] the lag-one correlation is close to zero\",\n      abs(lag_one_correlation(shuffled)) < 0.1)\nfor name, series in ((\"iid_lag_record.csv\", temp),\n                     (\"iid_lag_shuffled.csv\", shuffled)):\n    with open(OUT_DIR / name, \"w\") as f:\n        f.write(\"today,tomorrow\\n\")\n        for a, b in zip(series[:-1], series[1:]):\n            f.write(f\"{a:.1f},{b:.1f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-deseason]** Subtracting the seasonal average repairs the identically distributed half only: the lag-one correlation of what is left is still 0.67."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "WINDOW = 31\nwrapped = np.concatenate([temp[-(WINDOW // 2):], temp, temp[:WINDOW // 2]])\nseasonal = np.convolve(wrapped, np.ones(WINDOW) / WINDOW, mode=\"valid\")\nresidual = temp - seasonal\nprint(f\"[B-deseason] after subtracting the seasonal average the monthly \"\n      f\"means agree to {np.abs([residual[month == k].mean() for k in range(1, 13)]).max():.2f} deg, \"\n      f\"and the lag-one correlation is still \"\n      f\"{lag_one_correlation(residual):.3f}\")\ncheck(\"[B-deseason] the seasonal average is what separates January from July\",\n      abs(residual[month == 7].mean() - residual[month == 1].mean()) < 2.0)\ncheck(\"[B-deseason] the dependence survives it\",\n      lag_one_correlation(residual) > 0.5)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-verify]** The two standard methods applied to periods of growing length (January, January to June, the whole year): the largest gap between the empirical CDFs of two 15-day blocks of the period against the Kolmogorov-Smirnov threshold, and the lag-one correlation calibrated by a permutation test. Writes the panels of the entry's third figure."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "BLOCK = 15\nLEVEL = 0.05\nNR_PERM = 2000\nPERIODS = [(\"month\", \"January\", temp[month == 1], np.array(stamps)[month == 1]),\n           (\"halfyear\", \"January to June\", temp[month <= 6],\n            np.array(stamps)[month <= 6]),\n           (\"year\", \"the whole year\", temp, np.array(stamps))]\nverdict = {}\nfor tag, title, series, days in PERIODS:\n    nr_block = len(series) // BLOCK\n    blocks = [series[i * BLOCK:(i + 1) * BLOCK] for i in range(nr_block)]\n    spans = [f\"{days[i * BLOCK][5:10]} to {days[(i + 1) * BLOCK - 1][5:10]}\"\n             for i in range(nr_block)]\n    # the pair of blocks that are furthest apart, ties broken by the means\n    gap, spread, first, second = max(\n        (ks_statistic(blocks[i], blocks[j]),\n         abs(blocks[i].mean() - blocks[j].mean()), i, j)\n        for i in range(nr_block) for j in range(i + 1, nr_block))\n    corr = lag_one_correlation(series)\n    draws = np.array([abs(lag_one_correlation(rng.permutation(series)))\n                      for _ in range(NR_PERM)])\n    pvalue = (1 + (draws >= abs(corr)).sum()) / (NR_PERM + 1)\n    nr_pair = nr_block * (nr_block - 1) // 2\n    plain = ks_threshold(BLOCK, LEVEL)\n    adjusted = ks_threshold(BLOCK, LEVEL / nr_pair)\n    verdict[tag] = (title, len(series), nr_block, gap, corr, pvalue,\n                    (blocks[first], spans[first]),\n                    (blocks[second], spans[second]),\n                    nr_pair, plain, adjusted)\n    for label, index in ((\"early\", first), (\"late\", second)):\n        values, fraction = empirical_cdf(blocks[index])\n        with open(OUT_DIR / f\"iid_cdf_{tag}_{label}.csv\", \"w\") as f:\n            f.write(\"temp,frac\\n\")\n            for a, b in zip(values, fraction):\n                f.write(f\"{a:g},{b:.4f}\\n\")\n    print(f\"[B-verify] {title} ({len(series)} days, {nr_block} blocks of \"\n          f\"{BLOCK} days): largest gap between the empirical CDFs of two \"\n          f\"blocks {gap:.3f} ({spans[first]} against {spans[second]}); the \"\n          f\"correlation of consecutive values is {corr:.3f}, which no \"\n          f\"reordering among {NR_PERM} random ones reaches\")\n    print(f\"[B-verify] {title}: the Kolmogorov-Smirnov threshold at level \"\n          f\"{LEVEL} is {plain:.3f} for a single pair and {adjusted:.3f} after \"\n          f\"dividing the level among the {nr_pair} pairs; the gap \"\n          f\"{'exceeds' if gap > adjusted else 'stays below'} it\")\ncheck(\"[B-verify] the gap between two blocks grows with the length of the \"\n      \"period\", verdict[\"month\"][3] < verdict[\"halfyear\"][3] <= verdict[\"year\"][3])\ncheck(\"[B-verify] six months and a year separate two blocks completely\",\n      verdict[\"halfyear\"][3] == 1.0 and verdict[\"year\"][3] == 1.0)\ncheck(\"[B-verify] the lag-one correlation stays above 0.5 on every period\",\n      all(verdict[tag][4] > 0.5 for tag, _, _, _ in PERIODS))\ncheck(\"[B-verify] no permutation reaches the observed correlation\",\n      all(verdict[tag][5] < 2.0 / (NR_PERM + 1) for tag, _, _, _ in PERIODS))\ncheck(\"[B-verify] January stays below the threshold even for a single pair\",\n      verdict[\"month\"][3] < verdict[\"month\"][9])\ncheck(\"[B-verify] the longer periods exceed the threshold that accounts for \"\n      \"every pair of blocks\",\n      verdict[\"halfyear\"][3] > verdict[\"halfyear\"][10]\n      and verdict[\"year\"][3] > verdict[\"year\"][10])"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-optimal]** The same statistic on one month and on one year: the correlation of consecutive values over random reorderings, which tightens as the collection grows. Writes the panels of the entry's last figure."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "CASES = [(\"month\", \"January\", temp[month == 1]),\n         (\"year\", \"the whole year\", temp)]\noptimal = {}\nfor tag, title, series in CASES:\n    size = len(series)\n    corr = lag_one_correlation(series)\n    draws = np.array([lag_one_correlation(rng.permutation(series))\n                      for _ in range(NR_PERM)])\n    pvalue_corr = (1 + (draws >= corr).sum()) / (NR_PERM + 1)\n    optimal[tag] = (title, size, corr, draws, pvalue_corr)\n    counts, borders = np.histogram(draws, bins=28, range=(-0.65, 1.0))\n    with open(OUT_DIR / f\"iid_perm_{tag}.csv\", \"w\") as f:\n        f.write(\"center,frac\\n\")\n        for left, right, count in zip(borders[:-1], borders[1:], counts):\n            f.write(f\"{(left + right) / 2:.4f},{count / NR_PERM:.4f}\\n\")\n    print(f\"[B-optimal] {title} ({size} days): the correlation of \"\n          f\"consecutive values is {corr:.3f}, against random reorderings of \"\n          f\"spread {draws.std(ddof=1):.3f} \"\n          f\"(1/sqrt(m) is {1 / math.sqrt(size):.3f}), largest permuted value \"\n          f\"{draws.max():.3f}\")\ncheck(\"[B-optimal] the reordered correlations tighten as the period grows\",\n      optimal[\"year\"][3].std(ddof=1) < optimal[\"month\"][3].std(ddof=1))\ncheck(\"[B-optimal] their spread matches one over the square root of the \"\n      \"collection size\",\n      all(abs(optimal[tag][3].std(ddof=1) - 1 / math.sqrt(optimal[tag][1]))\n          < 0.03 for tag, _, _ in CASES))\ncheck(\"[B-optimal] no reordering reaches the observed correlation\",\n      all(optimal[tag][3].max() < optimal[tag][2] for tag, _, _ in CASES))\ncheck(\"[B-optimal] the reordered correlations stay far from the observed \"\n      \"value on both periods\",\n      all(optimal[tag][4] < 2.0 / (NR_PERM + 1) for tag, _, _ in CASES))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-plot]** preview"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "fig = plt.figure(figsize=(12, 12.5))\ngrid = fig.add_gridspec(4, 4, hspace=0.6, wspace=0.55)\nfor column, (_, tag, title) in enumerate(MONTHS):\n    ax = fig.add_subplot(grid[0, column])\n    day, values, _ = fine[tag]\n    ax.plot(day, values, color=\"0.3\", linewidth=0.5)\n    ax.set_xlabel(\"day of the month (UTC)\")\n    ax.set_ylabel(\"temperature in deg C\")\n    ax.set_ylim(-12, 38)\n    ax.set_title(f\"{title} 2024, every ten minutes\", fontsize=9)\nfor column, (series, title) in enumerate((\n        (temp, \"the record, in the order it was measured\"),\n        (shuffled, \"the same 366 numbers, permuted\"))):\n    ax = fig.add_subplot(grid[1, column])\n    ax.scatter(series[:-1], series[1:], s=9, color=\"0.45\")\n    ax.set_xlabel(\"daily maximum in deg C\")\n    ax.set_ylabel(\"next day in deg C\")\n    ax.set_title(title, fontsize=9)\n    ax.set_aspect(\"equal\")\nfor column, (tag, title, _, _) in enumerate(PERIODS):\n    ax = fig.add_subplot(grid[2, column])\n    for (block, span), style in ((verdict[tag][6], \"-\"),\n                                 (verdict[tag][7], \"--\")):\n        values, fraction = empirical_cdf(block)\n        ax.step(values, fraction, style, where=\"post\", color=\"0.3\",\n                linewidth=1.2, label=span)\n    ax.set_xlabel(\"daily maximum in deg C\")\n    ax.set_ylabel(\"fraction at or below\")\n    ax.set_title(f\"{title}: largest gap {verdict[tag][3]:.2f}\", fontsize=9)\n    ax.legend(frameon=False, fontsize=8)\nfor column, (tag, title, _) in enumerate(CASES):\n    ax = fig.add_subplot(grid[3, column])\n    ax.hist(optimal[tag][3], bins=28, range=(-0.65, 1.0), color=\"0.6\")\n    ax.axvline(optimal[tag][2], color=\"0.1\", linestyle=\"--\")\n    ax.annotate(\"observed\", (optimal[tag][2], 0), xytext=(-4, 12),\n                textcoords=\"offset points\", rotation=90, fontsize=8,\n                ha=\"right\")\n    ax.set_xlabel(\"correlation of consecutive values\")\n    ax.set_ylabel(\"reorderings\")\n    ax.set_title(f\"{title}: random reorderings\", fontsize=9)\nfig.suptitle(\"Krems 2024: the two halves of the i.i.d. property, checked on \"\n             \"months, on days, and on periods of growing length\", fontsize=11)\nfig.savefig(OUT_DIR / \"iid.png\", dpi=110, bbox_inches=\"tight\")\n\npassed = sum(ok for _, ok in report)\nprint(f\"\\n{passed}/{len(report)} checks pass\")\nif passed != len(report):\n    raise SystemExit(1)"
  }
 ]
}