{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "logreg.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# logistic regression \u2014 Python demo\n\nNumerical companion to the entry [logistic regression](https://dictionaryofml.org/terms/logreg.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nA one-feature binary classification trainset on which logistic regression is fit by a hand-written gradient descent loop: the average logistic loss decreases monotonically, the gradient vanishes at the learned parameters, the GD update leaves them (approximately) unchanged \u2014 its fixed point \u2014 and thresholding the learned hypothesis at zero classifies the trainset better than always answering with the majority label. Self-contained (numpy/matplotlib only), fixed seed.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/logreg.py`](https://dictionaryofml.org/terms/logreg.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(), \"logreg.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nlogreg.py \u2014 numerical companion to the glossary entry 'logistic regression'.\n\nPurpose\n-------\nA one-feature binary classification trainset on which logistic regression\nis fit by a hand-written gradient descent loop: the average logistic loss\ndecreases monotonically, the gradient vanishes at the learned parameters,\nthe GD update leaves them (approximately) unchanged \u2014 its fixed point \u2014\nand thresholding the learned hypothesis at zero classifies the trainset\nbetter than always answering with the majority label.  Self-contained\n(numpy/matplotlib only), fixed seed.\n\nBlocks\n------\n[B-data]    m = 40 data points with a single feature x drawn uniformly\n            from [-3, 3]; the binary label y in {-1, +1} is +1 with\n            probability sigmoid(2 x - 1), so the labels are noisy around\n            the point where 2 x - 1 = 0.\n[B-gd]      Gradient descent on the average logistic loss\n            f(w) = (1/m) sum_r log(1 + exp(-y^(r) w^T x^(r))) with the\n            constant feature 1 absorbing the intercept: the loss never\n            increases along the run, the gradient norm at the learned\n            parameters is small, and one further GD update moves them by\n            (approximately) nothing \u2014 the learned parameters are a fixed\n            point of the update.\n[B-clf]     The classifier sign(w^T x) obtained by thresholding the\n            learned hypothesis classifies a larger fraction of the\n            trainset correctly than the constant rule that always\n            answers with the majority label.\n[B-csv]     The two CSVs the entry's pgfplots figure reads.\n[B-preview] The matplotlib preview of the figure.\n\nOutputs\n-------\nlogreg_points.csv : the trainset, columns x,y01,cls (y01 = (y+1)/2 for\n                    drawing labels at heights 0 and 1; cls in {pos,neg}).\nlogreg_curve.csv  : the fitted probability curve sigmoid(w^T x) on a\n                    grid, columns x,p.\nlogreg.png        : matplotlib preview of the figure (checking only).\n\"\"\"\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}\")\n\n\ndef sigmoid(z):\n    return 1.0 / (1.0 + np.exp(-z))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-data]** m = 40 data points with a single feature x drawn uniformly from [-3, 3]; the binary label y in {-1, +1} is +1 with probability sigmoid(2 x - 1), so the labels are noisy around the point where 2 x - 1 = 0."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "rng = np.random.default_rng(3)\nm = 40\nx = rng.uniform(-3.0, 3.0, m)\np_true = sigmoid(2.0 * x - 1.0)\ny = np.where(rng.uniform(size=m) < p_true, 1.0, -1.0)\nX = np.c_[x, np.ones(m)]                    # constant feature absorbs the intercept\ncheck(\"[B-data]    both labels occur\", (y > 0).any() and (y < 0).any())"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-gd]** Gradient descent on the average logistic loss f(w) = (1/m) sum_r log(1 + exp(-y^(r) w^T x^(r))) with the constant feature 1 absorbing the intercept: the loss never increases along the run, the gradient norm at the learned parameters is small, and one further GD update moves them by (approximately) nothing \u2014 the learned parameters are a fixed point of the update."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def avg_logloss(w):\n    return float(np.mean(np.log1p(np.exp(-y * (X @ w)))))\n\n\ndef grad(w):\n    return -(X * (y * sigmoid(-y * (X @ w)))[:, None]).mean(axis=0)\n\n\neta = 0.5\nw = np.zeros(2)\nlosses = [avg_logloss(w)]\nfor _ in range(2000):\n    w = w - eta * grad(w)                   # the GD update w <- T(w)\n    losses.append(avg_logloss(w))\ncheck(\"[B-gd]      the average logistic loss never increases and ends \"\n      \"below its start\",\n      all(b <= a + 1e-12 for a, b in zip(losses, losses[1:]))\n      and losses[-1] < losses[0])\ncheck(\"[B-gd]      the gradient vanishes at the learned parameters\",\n      float(np.linalg.norm(grad(w))) < 1e-4)\nw_next = w - eta * grad(w)                  # one further update\ncheck(\"[B-gd]      the learned parameters are a fixed point of the update\",\n      float(np.linalg.norm(w_next - w)) < 1e-4)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-clf]** The classifier sign(w^T x) obtained by thresholding the learned hypothesis classifies a larger fraction of the trainset correctly than the constant rule that always answers with the majority label."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "frac_lr = float(np.mean(np.sign(X @ w) == y))\nmajority = 1.0 if (y > 0).sum() >= (y < 0).sum() else -1.0\nfrac_const = float(np.mean(y == majority))\nprint(f\"    correctly classified: logistic regression {frac_lr:.2f}, \"\n      f\"majority label {frac_const:.2f}\")\ncheck(\"[B-clf]     sign(w^T x) beats the constant majority rule\",\n      frac_lr > frac_const)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-csv]** The two CSVs the entry's pgfplots figure reads."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "with open(OUT_DIR / \"logreg_points.csv\", \"w\") as f:\n    f.write(\"x,y01,cls\\n\")\n    for xi, yi in zip(x, y):\n        f.write(f\"{xi:.4f},{(yi + 1) / 2:.0f},{'pos' if yi > 0 else 'neg'}\\n\")\n\ngx = np.linspace(-3.2, 3.2, 161)\ngp = sigmoid(w[0] * gx + w[1])\nwith open(OUT_DIR / \"logreg_curve.csv\", \"w\") as f:\n    f.write(\"x,p\\n\")\n    for xi, pi in zip(gx, gp):\n        f.write(f\"{xi:.4f},{pi:.4f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-preview]** The matplotlib preview of the figure."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "fig, ax = plt.subplots(figsize=(5.2, 3.4))\nax.plot(gx, gp, \"k-\", lw=1.6, label=\"sigmoid($\\\\hat{w}^{\\\\top} x$)\")\nax.plot(x[y > 0], np.ones((y > 0).sum()), \"ko\", ms=5, label=\"label $+1$\")\nax.plot(x[y < 0], np.zeros((y < 0).sum()), \"ks\", mfc=\"none\", ms=5,\n        label=\"label $-1$\")\nax.set_xlabel(\"feature $x$\")\nax.set_ylabel(\"label / probability of label $+1$\")\nax.set_title(\"logistic regression: fitted probability curve\")\nax.legend(frameon=False, loc=\"center right\", fontsize=8)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"logreg.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 / 'logreg_points.csv'}, {OUT_DIR / 'logreg_curve.csv'}, \"\n      f\"{OUT_DIR / 'logreg.png'}\")\nif n_ok != len(report):\n    raise SystemExit(1)"
  }
 ]
}