{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "samplesize.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# sample size \u2014 Python demo\n\nNumerical companion to the entry [sample size](https://dictionaryofml.org/terms/samplesize.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nBacks the entry's chain of claims on the frost-warning hypothesis of Krems an der Donau: a hypothesis is evaluated by its expected loss (the risk) under the i.i.d. assumption, a method sees only a sample and estimates the risk by a sample average, the gap between the two shrinks with the sample size at the rate Hoeffding's inequality states, and that inequality turns a required accuracy and confidence into a sample size. Self- contained (numpy/matplotlib only, stdlib urllib for the download), fixed seed.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/samplesize.py`](https://dictionaryofml.org/terms/samplesize.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(), \"samplesize.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nsamplesize.py \u2014 numerical companion to the glossary entry 'sample size'.\n\nPurpose\n-------\nBacks the entry's chain of claims on the frost-warning hypothesis of Krems\nan der Donau: a hypothesis is evaluated by its expected loss (the risk)\nunder the i.i.d. assumption, a method sees only a sample and estimates\nthe risk by a sample average, the gap between the two shrinks with the\nsample size at the rate Hoeffding's inequality states, and that inequality\nturns a required accuracy and confidence into a sample size.  Self-\ncontained (numpy/matplotlib only, stdlib urllib for the download), fixed\nseed.\n\nSetup\n-----\nDaily minimum and maximum temperature of 2024 at Krems (GeoSphere\nAustria).  A data point is one day: features (minimum, maximum), label 1\nif the next day's minimum is below 0 degrees (frost) and -1 otherwise.\nThe hypothesis under evaluation is fixed before any sampling: it predicts\nfrost whenever the day's minimum is below 1 degree.  The 0/1 loss of this\nhypothesis on the 365 days defines the population; its average over all\ndays plays the role of the risk, and a sample of size m is m days drawn\nuniformly with replacement, so the sampled losses are i.i.d. RVs.\n\nBlocks\n------\n[B-risk]    The hypothesis, its 0/1 loss on every day, and the risk (the\n            loss averaged over the whole year).\n[B-sample]  Samples of size m = 10, 30, 100, 300: the sample average of\n            the loss estimates the risk; over 2000 draws the largest\n            deviation shrinks as m grows, and the spread of the estimate\n            falls like 1/sqrt(m).\n[B-bound]   Hoeffding: the fraction of draws with a deviation above\n            epsilon never exceeds 2 exp(-2 m epsilon^2); the bound is\n            loose but has the right shape.\n[B-choose]  Sample size from the bound: a deviation of at most 0.05 with\n            probability at least 0.95 needs m >= log(2/0.05) / (2 0.05^2)\n            = 738 days, and the empirical rate at m = 738 confirms it.\n[B-dependent] Consecutive days are not i.i.d.: the 0/1 loss of the\n            hypothesis has sample autocorrelation 0.24 at lag 1 (a wrong\n            prediction follows a wrong prediction).  The variance of the\n            average over m consecutive days is that over m random days\n            times 1 + 2 sum_k (1 - k/m) rho_k; for m = 30 the factor is\n            4.2 and the measured ratio 4.4, so 30 consecutive days are\n            worth 7 random ones: the effective sample size.\n[B-corr]    Sample size for a statistical test: are two features of a day,\n            its minimum temperature and its diurnal range (maximum minus\n            minimum), correlated?  Over the year the Pearson correlation\n            coefficient is 0.16.  With the features uncorrelated, the\n            sample correlation r scatters around 0 with standard deviation\n            1/sqrt(m) and the test \"reject 'uncorrelated' when\n            |artanh(r)| sqrt(m - 3) > 1.96\" rejects with frequency 0.05;\n            detecting 0.16 with probability 0.8 needs\n            m >= 3 + ((1.96 + 0.84) / artanh(0.16))^2 = 304 days (Fisher's\n            formula for Gaussian features): the empirical power is 0.84\n            there and 0.10 at m = 30.\n\nOutputs\n-------\nsamplesize_losses.csv   : day, minimum, maximum, label, loss of the hypothesis\nsamplesize_curve.csv    : m, sample standard deviation of the estimate,\n                          1/(2 sqrt m), Hoeffding width for delta = 0.05\nsamplesize_bound.csv    : m, epsilon, empirical tail probability, Hoeffding\n                          bound\nsamplesize_corr.csv     : m, standard deviation of artanh(r), 1/sqrt(m - 3),\n                          power of the correlation test at level 0.05\nsamplesize_acf.csv      : lag, sample autocorrelation of the loss sequence\nsamplesize.png          : matplotlib preview (checking only)\n\"\"\"\n\nimport json\nimport math\nimport urllib.request\n\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nfrom pathlib import Path\n\nOUT_DIR = Path(__file__).parent\n\nreport = []\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-risk]** The hypothesis, its 0/1 loss on every day, and the risk (the loss averaged over the whole year)."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "URL = (\"https://dataset.api.hub.geosphere.at/v1/station/historical/\"\n       \"klima-v2-1d?parameters=tlmin,tlmax&station_ids=3805\"\n       \"&start=2024-01-01&end=2024-12-31\")\nwith urllib.request.urlopen(URL, timeout=120) as resp:\n    payload = json.load(resp)\npar = payload[\"features\"][0][\"properties\"][\"parameters\"]\ntmin = np.array(par[\"tlmin\"][\"data\"], dtype=float)\ntmax = np.array(par[\"tlmax\"][\"data\"], dtype=float)\ndays = [s[:10] for s in payload[\"timestamps\"]]\nx = np.c_[tmin[:-1], tmax[:-1]]                 # a day's two features\ny = np.where(tmin[1:] < 0.0, 1, -1)             # frost on the following day\nN = len(y)\nTHRESH = 1.0\npred = np.where(x[:, 0] < THRESH, 1, -1)        # the fixed hypothesis\nloss = (pred != y).astype(float)                # 0/1 loss per day\nrisk = float(loss.mean())\ncheck(f\"[B-risk]    {N} days of 2024, {int((y == 1).sum())} followed by frost; the hypothesis \"\n      f\"'frost if today's minimum is below {THRESH:g} degrees' has risk {risk:.3f} \"\n      f\"(its 0/1 loss averaged over the year)\", 0.0 < risk < 0.5)\nwith open(OUT_DIR / \"samplesize_losses.csv\", \"w\") as fh:\n    fh.write(\"day,tmin,tmax,label,loss\\n\")\n    for d, xi, yi, li in zip(days[:-1], x, y, loss):\n        fh.write(f\"{d},{xi[0]:g},{xi[1]:g},{yi},{int(li)}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-sample]** Samples of size m = 10, 30, 100, 300: the sample average of the loss estimates the risk; over 2000 draws the largest deviation shrinks as m grows, and the spread of the estimate falls like 1/sqrt(m)."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "rng = np.random.default_rng(2024)\nSIZES = [10, 30, 100, 300, 1000]\nDRAWS = 2000\nest = {m: loss[rng.integers(0, N, (DRAWS, m))].mean(axis=1) for m in SIZES}\ndev = {m: np.abs(est[m] - risk) for m in SIZES}\nsd = {m: float(est[m].std()) for m in SIZES}\nprint(\"    sample size m, largest |estimate - risk| over 2000 draws, spread of the estimate:\")\nfor m in SIZES:\n    print(f\"      m = {m:5d}: largest deviation {dev[m].max():.3f}, standard deviation {sd[m]:.4f}\")\ncheck(\"[B-sample]  the sample average of the loss estimates the risk, and its largest \"\n      \"deviation over 2000 draws shrinks with m\",\n      all(dev[a].max() > dev[b].max() for a, b in zip(SIZES, SIZES[1:])))\nratio = [sd[m] * math.sqrt(m) for m in SIZES]\ncheck(f\"[B-sample]  the spread of the estimate falls like 1/sqrt(m): sd * sqrt(m) stays near \"\n      f\"{np.mean(ratio):.2f} ({min(ratio):.2f} to {max(ratio):.2f})\",\n      max(ratio) - min(ratio) < 0.1)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-bound]** Hoeffding: the fraction of draws with a deviation above epsilon never exceeds 2 exp(-2 m epsilon^2); the bound is loose but has the right shape."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "EPS = [0.05, 0.1, 0.15]\nrows = []\nfor m in SIZES:\n    for eps in EPS:\n        emp = float(np.mean(dev[m] > eps))\n        bnd = min(1.0, 2.0 * math.exp(-2.0 * m * eps ** 2))\n        rows.append((m, eps, emp, bnd))\ncheck(\"[B-bound]   Hoeffding: the fraction of draws with a deviation above epsilon never exceeds \"\n      \"2 exp(-2 m epsilon^2), for every m and epsilon\", all(emp <= bnd + 1e-12 for _, _, emp, bnd in rows))\nwith open(OUT_DIR / \"samplesize_bound.csv\", \"w\") as fh:\n    fh.write(\"m,epsilon,empirical,hoeffding\\n\")\n    for m, eps, emp, bnd in rows:\n        fh.write(f\"{m},{eps},{emp:.4f},{bnd:.4f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-choose]** Sample size from the bound: a deviation of at most 0.05 with probability at least 0.95 needs m >= log(2/0.05) / (2 0.05^2) = 738 days, and the empirical rate at m = 738 confirms it."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "EPS0, DELTA = 0.05, 0.05\nm_needed = math.ceil(math.log(2.0 / DELTA) / (2.0 * EPS0 ** 2))\nest_needed = loss[rng.integers(0, N, (DRAWS, m_needed))].mean(axis=1)\nrate = float(np.mean(np.abs(est_needed - risk) <= EPS0))\ncheck(f\"[B-choose]  a deviation of at most {EPS0} with probability at least {1 - DELTA} needs \"\n      f\"m >= log(2/{DELTA}) / (2 * {EPS0}^2) = {m_needed} days; at that size \"\n      f\"{rate:.3f} of 2000 draws are within {EPS0}\", rate >= 1 - DELTA)\ngrid = np.unique(np.round(np.logspace(1, 3.2, 40)).astype(int))\nwith open(OUT_DIR / \"samplesize_curve.csv\", \"w\") as fh:\n    fh.write(\"m,sd,halfsqrt,hoeffding\\n\")\n    for m in grid:\n        s = float(loss[rng.integers(0, N, (400, m))].mean(axis=1).std())\n        fh.write(f\"{m},{s:.4f},{1 / (2 * math.sqrt(m)):.4f},{math.sqrt(math.log(2 / DELTA) / (2 * m)):.4f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-dependent]** Consecutive days are not i.i.d.: the 0/1 loss of the hypothesis has sample autocorrelation 0.24 at lag 1 (a wrong prediction follows a wrong prediction). The variance of the average over m consecutive days is that over m random days times 1 + 2 sum_k (1 - k/m) rho_k; for m = 30 the factor is 4.2 and the measured ratio 4.4, so 30 consecutive days are worth 7 random ones: the effective sample size."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def acf(v, k):\n    v = v - v.mean()\n    return float((v[:-k] * v[k:]).sum() / (v * v).sum())\n\n\nrhos = [acf(loss, k) for k in range(1, 61)]\nprint(\"    sample autocorrelation of the loss at lags 1 to 7: \"\n      + \", \".join(f\"{r:.2f}\" for r in rhos[:7]))\ndep = {}\nfor m in (30, 60):\n    windows = np.array([loss[i:i + m].mean() for i in range(N - m + 1)])   # consecutive days\n    random_ = loss[rng.integers(0, N, (4000, m))].mean(1)                   # days drawn at random\n    factor = 1 + 2 * sum((1 - k / m) * rhos[k - 1] for k in range(1, m))\n    dep[m] = (windows.var() / random_.var(), factor, m / factor)\n    print(f\"      m = {m}: variance ratio consecutive / random {dep[m][0]:.2f}, \"\n          f\"factor 1 + 2 sum (1 - k/m) rho_k = {factor:.2f}, effective sample size {m / factor:.0f}\")\ncheck(f\"[B-dependent] the loss of consecutive days is autocorrelated, rho_1 = {rhos[0]:.2f}, \"\n      f\"and the average over m consecutive days has the variance of m random days times \"\n      f\"1 + 2 sum (1 - k/m) rho_k, to within 25% at m = 30 and 60\",\n      rhos[0] > 0.15 and all(abs(r / f - 1) < 0.25 for r, f, _ in dep.values()))\ncheck(f\"[B-dependent] 30 consecutive days are worth {dep[30][2]:.0f} random ones (effective sample size)\",\n      5 <= dep[30][2] <= 9)\nwith open(OUT_DIR / \"samplesize_acf.csv\", \"w\") as fh:\n    fh.write(\"lag,rho\\n\")\n    for k, r in enumerate(rhos[:30], 1):\n        fh.write(f\"{k},{r:.4f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-corr]** Sample size for a statistical test: are two features of a day, its minimum temperature and its diurnal range (maximum minus minimum), correlated? Over the year the Pearson correlation coefficient is 0.16. With the features uncorrelated, the sample correlation r scatters around 0 with standard deviation 1/sqrt(m) and the test \"reject 'uncorrelated' when |artanh(r)| sqrt(m - 3) > 1.96\" rejects with frequency 0.05; detecting 0.16 with probability 0.8 needs m >= 3 + ((1.96 + 0.84) / artanh(0.16))^2 = 304 days (Fisher's formula for Gaussian features): the empirical power is 0.84 there and 0.10 at m = 30."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def gauss_quantile(p):\n    \"\"\"The standard Gaussian quantile, by bisection on the cdf.\"\"\"\n    lo, hi = -10.0, 10.0\n    for _ in range(100):\n        mid = 0.5 * (lo + hi)\n        if 0.5 * (1.0 + math.erf(mid / math.sqrt(2.0))) < p:\n            lo = mid\n        else:\n            hi = mid\n    return 0.5 * (lo + hi)\n\n\nALPHA, POWER = 0.05, 0.8\nz_alpha, z_power = gauss_quantile(1 - ALPHA / 2), gauss_quantile(POWER)\nf1, f2 = x[:, 0], x[:, 1] - x[:, 0]                # two features of the 365 labeled days\nrho = float(np.corrcoef(f1, f2)[0, 1])\nm_corr = math.ceil(3 + ((z_alpha + z_power) / math.atanh(rho)) ** 2)\n\n\ndef sample_corr(a, b, m, draws):\n    idx = rng.integers(0, len(a), (draws, m))\n    aa, bb = a[idx], b[idx]\n    aa = aa - aa.mean(1, keepdims=True); bb = bb - bb.mean(1, keepdims=True)\n    return (aa * bb).sum(1) / np.sqrt((aa ** 2).sum(1) * (bb ** 2).sum(1))\n\n\ndef power_at(m, draws=2000):\n    r = sample_corr(f1, f2, m, draws)\n    return float(np.mean(np.abs(np.arctanh(r)) * math.sqrt(m - 3) > z_alpha))\n\n\nprint(f\"    minimum temperature and diurnal range: correlation {rho:.3f} over the {len(f1)} days\")\n\n\ndef null_corr(m, draws):\n    \"\"\"Sample correlation of the two features drawn independently: uncorrelated by construction.\"\"\"\n    a = f1[rng.integers(0, len(f1), (draws, m))]; b = f2[rng.integers(0, len(f2), (draws, m))]\n    a = a - a.mean(1, keepdims=True); b = b - b.mean(1, keepdims=True)\n    return (a * b).sum(1) / np.sqrt((a ** 2).sum(1) * (b ** 2).sum(1))\n\n\nnsd = {m: float(null_corr(m, 2000).std() * math.sqrt(m)) for m in (30, 100, 300, 1000)}\ncheck(\"[B-corr]    with the two features uncorrelated, the sample correlation r scatters around 0 \"\n      \"with standard deviation close to 1/sqrt(m): sd * sqrt(m) = \"\n      + \", \".join(f\"{v:.2f} (m = {m})\" for m, v in nsd.items()), all(0.88 < v < 1.12 for v in nsd.values()))\nsize = float(np.mean(np.abs(np.arctanh(null_corr(100, 4000))) * math.sqrt(97) > z_alpha))\ncheck(f\"[B-corr]    on uncorrelated features the test rejects at level {ALPHA} with frequency \"\n      f\"{size:.3f}\", size < 0.075)\nzsd = {m: float(np.arctanh(sample_corr(f1, f2, m, 2000)).std() * math.sqrt(m - 3)) for m in (30, 300, 1000)}\nprint(\"    at rho = 0.16, sd of artanh(r) times sqrt(m - 3): \"\n      + \", \".join(f\"{v:.2f} (m = {m})\" for m, v in zsd.items()) + \" (1 for Gaussian features)\")\np_small, p_need = power_at(30), power_at(m_corr)\ncheck(f\"[B-corr]    detecting rho = {rho:.2f} at level {ALPHA} with power {POWER} needs \"\n      f\"m >= 3 + (({z_alpha:.2f} + {z_power:.2f}) / artanh({rho:.2f}))^2 = {m_corr} days; empirical power \"\n      f\"{p_need:.2f} at m = {m_corr}, only {p_small:.2f} at m = 30\",\n      abs(p_need - POWER) < 0.08 and p_small < 0.3)\nwith open(OUT_DIR / \"samplesize_corr.csv\", \"w\") as fh:\n    fh.write(\"m,zsd,invsqrt,power\\n\")\n    for m in np.unique(np.round(np.logspace(1, 3.2, 30)).astype(int)):\n        r = sample_corr(f1, f2, int(m), 1000)\n        fh.write(f\"{m},{np.arctanh(r).std():.4f},{1 / math.sqrt(m - 3):.4f},\"\n                 f\"{np.mean(np.abs(np.arctanh(r)) * math.sqrt(m - 3) > z_alpha):.3f}\\n\")\n\n# -------------------------------------------------------------- preview\nfig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(13.5, 3.7))\nfor m, c in zip(SIZES[:4], (\"0.85\", \"0.6\", \"0.35\", \"0.0\")):\n    ax1.hist(est[m], bins=np.linspace(0, 0.5, 51), histtype=\"step\", color=c, label=f\"m = {m}\")\nax1.axvline(risk, color=\"k\", ls=\"--\", lw=0.8, label=f\"risk {risk:.3f}\")\nax1.set_xlabel(\"sample average of the 0/1 loss\"); ax1.set_ylabel(\"draws out of 2000\")\nax1.set_title(\"the estimate concentrates as m grows\"); ax1.legend(frameon=False, fontsize=8)\ncv = np.genfromtxt(OUT_DIR / \"samplesize_curve.csv\", delimiter=\",\", names=True)\nax2.loglog(cv[\"m\"], cv[\"sd\"], \"k-\", label=\"spread of the estimate\")\nax2.loglog(cv[\"m\"], cv[\"halfsqrt\"], \"k:\", label=r\"$1/(2\\sqrt{m})$\")\nax2.loglog(cv[\"m\"], cv[\"hoeffding\"], \"k--\", label=r\"Hoeffding width, $\\delta = 0.05$\")\nax2.axvline(m_needed, color=\"0.6\", lw=0.8); ax2.text(m_needed * 1.05, 0.3, f\"m = {m_needed}\", fontsize=8)\nax2.set_xlabel(\"sample size $m$\"); ax2.set_ylabel(\"deviation from the risk\")\nax2.set_title(\"choosing the sample size\"); ax2.legend(frameon=False, fontsize=8)\ncc = np.genfromtxt(OUT_DIR / \"samplesize_corr.csv\", delimiter=\",\", names=True)\nax3.semilogx(cc[\"m\"], cc[\"power\"], \"k-\", label=f\"power of the test, correlation {rho:.2f}\")\nax3.axhline(POWER, color=\"0.6\", lw=0.8, ls=\"--\", label=f\"power {POWER}\")\nax3.axvline(m_corr, color=\"0.6\", lw=0.8); ax3.text(m_corr * 1.08, 0.1, f\"m = {m_corr}\", fontsize=8)\nax3.set_xlabel(\"sample size $m$\"); ax3.set_ylabel(\"probability of rejecting 'uncorrelated'\")\nax3.set_title(\"sample size for a correlation test\"); ax3.legend(frameon=False, fontsize=8)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"samplesize.png\", dpi=110)\n\nn_ok = sum(ok for _, ok in report)\nprint(f\"\\n{n_ok}/{len(report)} checks pass\")\nprint(f\"wrote {OUT_DIR / 'samplesize_losses.csv'}, {OUT_DIR / 'samplesize_curve.csv'}, \"\n      f\"{OUT_DIR / 'samplesize_bound.csv'}, {OUT_DIR / 'samplesize_corr.csv'}, \"\n      f\"{OUT_DIR / 'samplesize_acf.csv'}, {OUT_DIR / 'samplesize.png'}\")\nif n_ok != len(report):\n    raise SystemExit(1)"
  }
 ]
}