{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "gaussrv.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# Gaussian random variable (Gaussian RV) \u2014 Python demo\n\nNumerical companion to the entry [Gaussian random variable (Gaussian RV)](https://dictionaryofml.org/terms/gaussrv.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nOne block per paragraph of the entry (marked [P...]), in entry order: each block verifies numerically what the corresponding paragraph asserts. Self-contained (numpy/matplotlib only, math.erf for the CDF), fixed seed.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/gaussrv.py`](https://dictionaryofml.org/terms/gaussrv.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(), \"gaussrv.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\ngaussrv.py \u2014 numerical companion to the glossary entry 'gaussrv'.\n\nOne block per paragraph of the entry (marked [P...]), in entry order: each\nblock verifies numerically what the corresponding paragraph asserts.\nSelf-contained (numpy/matplotlib only, math.erf for the CDF), fixed seed.\n\nBlocks\n------\n[P-def]   A standard Gaussian RV z has the pdf p(eta) = exp(-eta^2/2)/\n          sqrt(2 pi): it integrates to 1 and peaks at 1/sqrt(2 pi). The\n          general Gaussian RV x := sigma z + mu has mean mu and variance\n          sigma^2, and the fraction of its realizations below eta matches\n          Phi((eta - mu)/sigma).\n[P-gen]   The two constructions of the entry. Exact: mu + sigma\n          Phi^{-1}(U) with U uniform on (0,1) reproduces N(mu, sigma^2).\n          Approximate: the normalized sum of n i.i.d. coin flips, each\n          +1 or -1 with probability 1/2 (mean 0, variance 1). Its exact\n          probabilities are computed from the number of head patterns,\n          so the largest gap to the Gaussian CDF is exact, not measured:\n          it shrinks from 0.34 at n = 1 to 0.05 at n = 64 - and stays\n          positive, because a finite sum takes only n + 1 values and is\n          therefore not Gaussian.\n[P-param] mu and sigma act separately on the pdf: adding mu shifts it\n          along the horizontal axis without changing its shape, while\n          multiplying by sigma > 0 stretches it by the factor sigma and\n          divides its peak by sigma, so the area stays 1. The mean and\n          the variance read off the pdf by integration are mu and\n          sigma^2, and the width at half the peak is proportional to\n          sigma.\n[P-vec]   A Gaussian random vector x := A z + mu, with z a vector of\n          i.i.d. standard Gaussian RVs and A a matrix square root of C\n          (a Cholesky factor), has mean mu and covariance matrix C. Its\n          components are independent exactly when C is diagonal: for a\n          diagonal C the joint probabilities factorize into the\n          marginals, for a non-diagonal one they do not.\n[P-gp]    A Gaussian random vector is a stochastic process indexed by\n          {1, ..., d}: its restriction to a subset of the indices is\n          again a Gaussian random vector, with the covariance matrix cut\n          out of C along those indices.\n[P-conc]  A Lipschitz function of a standard Gaussian random vector\n          concentrates: the Euclidean norm has Lipschitz constant 1, so\n          the length of z in d dimensions has expectation at most\n          sqrt(d), variance at most 1, and leaves a window of width t\n          around its expectation with probability at most\n          2 exp(-t^2/2) \u2014 at d = 10, 100 and 1000 alike. The squared\n          length, which is not Lipschitz, instead fluctuates ten times\n          more at d = 1000 than at d = 10. The bound constrains the\n          tails only: the largest of 50 standard Gaussian RVs obeys it\n          and is still skewed, hence not Gaussian.\n[P-ent]   Among RVs with a given variance, the Gaussian one maximizes\n          the differential entropy: by numerical integration its entropy\n          is log(2 pi e sigma^2)/2, above that of a uniform RV and of a\n          sum of two uniform RVs with the same variance.\n[P-reg]   For the regression model y = w^T x + eps with eps a Gaussian\n          RV of variance sigma^2, the negative log-likelihood of the\n          labels of the training set is an increasing affine function of\n          the average squared error loss. Maximizing the likelihood over\n          w therefore returns the same weights as minimizing that loss,\n          which is the ERM problem of linear regression.\n\nOutputs\n-------\ngaussrv_clt.csv    : the exact probability density of the normalized sum\n                     of 16 coin flips, next to the standard Gaussian pdf\n                     at the same points (s, density, gausspdf).\ngaussrv_cltgap.csv : the largest gap between the CDF of that sum and the\n                     Gaussian CDF, against the number of coin flips\n                     (n, gap).\ngaussrv.png        : preview figure (checking only).\n\nData generated by pythondemos/gaussrv.py.\n\"\"\"\n\nimport math\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\nfrom pathlib import Path\n\nOUT_DIR = Path(__file__).parent\n\nrng = np.random.default_rng(42)\nreport = []\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")\n\n\n_ERF = np.frompyfunc(math.erf, 1, 1)\n\n\ndef Phi(t):\n    \"\"\"CDF of a standard Gaussian RV.\"\"\"\n    t = np.asarray(t, dtype=float)\n    return 0.5 * (1.0 + np.asarray(_ERF(t / np.sqrt(2.0)), dtype=float))\n\n\n_T_GRID = np.linspace(-8.0, 8.0, 200_001)\n_P_GRID = Phi(_T_GRID)\n\n\ndef Phi_inv(u):\n    \"\"\"Inverse of Phi, by interpolation of Phi on a fine grid.\"\"\"\n    return np.interp(u, _P_GRID, _T_GRID)\n\n\ndef pdf_std(eta):\n    \"\"\"pdf of a standard Gaussian RV.\"\"\"\n    eta = np.asarray(eta, dtype=float)\n    return np.exp(-(eta**2) / 2.0) / np.sqrt(2.0 * np.pi)\n\n\ndef pdf_gen(eta, mu, sigma):\n    \"\"\"pdf of a Gaussian RV with mean mu and variance sigma^2.\"\"\"\n    return pdf_std((np.asarray(eta, dtype=float) - mu) / sigma) / sigma"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-def]** A standard Gaussian RV z has the pdf p(eta) = exp(-eta^2/2)/ sqrt(2 pi): it integrates to 1 and peaks at 1/sqrt(2 pi). The general Gaussian RV x := sigma z + mu has mean mu and variance sigma^2, and the fraction of its realizations below eta matches Phi((eta - mu)/sigma)."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-def] the standard Gaussian pdf, and x := sigma z + mu\")\ngrid = np.linspace(-12.0, 12.0, 240_001)\ncheck(\"the standard Gaussian pdf integrates to 1\",\n      abs(np.trapezoid(pdf_std(grid), grid) - 1.0) < 1e-9)\ncheck(\"its peak is 1/sqrt(2 pi) = 0.3989\",\n      abs(pdf_std(0.0) - 1.0 / np.sqrt(2.0 * np.pi)) < 1e-12)\nmu, sigma = 2.0, 1.5\nm = 10**6\nz = rng.standard_normal(m)\nx = sigma * z + mu\ncheck(\"sample mean of x := sigma z + mu matches mu = 2\",\n      abs(x.mean() - mu) < 1e-2)\ncheck(\"sample variance of x matches sigma^2 = 2.25\",\n      abs(np.mean((x - x.mean()) ** 2) - sigma**2) < 2e-2)\nfor eta in (-1.0, 2.0, 4.5):\n    frac = np.mean(x <= eta)\n    check(f\"fraction of realizations below {eta} matches \"\n          f\"Phi((eta - mu)/sigma)\",\n          abs(frac - float(Phi((eta - mu) / sigma))) < 3e-3)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-gen]** The two constructions of the entry. Exact: mu + sigma Phi^{-1}(U) with U uniform on (0,1) reproduces N(mu, sigma^2). Approximate: the normalized sum of n i.i.d. coin flips, each +1 or -1 with probability 1/2 (mean 0, variance 1). Its exact probabilities are computed from the number of head patterns, so the largest gap to the Gaussian CDF is exact, not measured: it shrinks from 0.34 at n = 1 to 0.05 at n = 64 - and stays positive, because a finite sum takes only n + 1 values and is therefore not Gaussian."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-gen] exact generation by the inverse CDF, approximate \"\n      \"generation by summing coin flips\")\nu = rng.random(200_000)\ncheck(\"Phi(Phi^{-1}(u)) returns u\", np.max(np.abs(Phi(Phi_inv(u)) - u)) < 1e-7)\nx_inv = mu + sigma * Phi_inv(u)\ncheck(\"sample mean of mu + sigma Phi^{-1}(U) matches mu\",\n      abs(x_inv.mean() - mu) < 2e-2)\ncheck(\"its sample variance matches sigma^2\",\n      abs(np.mean((x_inv - x_inv.mean()) ** 2) - sigma**2) < 5e-2)\ncheck(\"the fraction of its realizations below 2.5 matches \"\n      \"Phi((2.5 - mu)/sigma)\",\n      abs(np.mean(x_inv <= 2.5) - float(Phi((2.5 - mu) / sigma))) < 5e-3)\n\n\ndef coinflip_sum(n):\n    \"\"\"Exact probabilities of the normalized sum of n coin flips.\n\n    Each flip is +1 or -1 with probability 1/2, so the sum has mean 0 and\n    variance n; k heads give the sum 2k - n. The number of flip patterns\n    with k heads is the count of k-subsets of the n flips.\n    \"\"\"\n    k = np.arange(n + 1)\n    patterns = np.array([math.comb(n, int(kk)) for kk in k], dtype=float)\n    prob = patterns / 2.0**n\n    s = (2.0 * k - n) / np.sqrt(n)\n    return s, prob\n\n\ndef cdf_gap(n):\n    \"\"\"Largest gap between the CDF of the normalized sum and Phi.\"\"\"\n    s, prob = coinflip_sum(n)\n    upper = np.cumsum(prob)               # CDF at and above each value\n    lower = upper - prob                  # CDF just below each value\n    target = Phi(s)\n    return float(max(np.max(np.abs(upper - target)),\n                     np.max(np.abs(lower - target))))\n\n\ns16, p16 = coinflip_sum(16)\ncheck(\"the normalized sum of 16 coin flips has mean 0\",\n      abs(float(np.sum(p16 * s16))) < 1e-12)\ncheck(\"it has variance 1\", abs(float(np.sum(p16 * s16**2)) - 1.0) < 1e-12)\ncheck(\"it takes 17 values, so it is not Gaussian\", s16.size == 16 + 1)\ngaps = {n: cdf_gap(n) for n in (1, 4, 16, 64)}\ncheck(\"the gap to the Gaussian CDF shrinks with n\",\n      gaps[1] > gaps[4] > gaps[16] > gaps[64])\ncheck(\"at n = 1 the sum is a single coin flip and the gap is \"\n      \"Phi(1) - 1/2 = 0.3413\",\n      abs(gaps[1] - (float(Phi(1.0)) - 0.5)) < 1e-12)\ncheck(\"at n = 64 the gap is below 0.05\", gaps[64] < 0.05)\ncheck(\"every finite n leaves a positive gap\",\n      all(g > 0.0 for g in gaps.values()))\nn_grid = np.arange(1, 65)\ngap_grid = np.array([cdf_gap(int(n)) for n in n_grid])\nnp.savetxt(OUT_DIR / \"gaussrv_cltgap.csv\",\n           np.column_stack([n_grid, gap_grid]),\n           fmt=[\"%d\", \"%.6f\"], delimiter=\",\", header=\"n,gap\", comments=\"\")\nwidth16 = 2.0 / np.sqrt(16)               # spacing of the 17 values\nnp.savetxt(OUT_DIR / \"gaussrv_clt.csv\",\n           np.column_stack([s16, p16 / width16, pdf_std(s16)]),\n           fmt=[\"%.4f\", \"%.6f\", \"%.6f\"], delimiter=\",\",\n           header=\"s,density,gausspdf\", comments=\"\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-param]** mu and sigma act separately on the pdf: adding mu shifts it along the horizontal axis without changing its shape, while multiplying by sigma > 0 stretches it by the factor sigma and divides its peak by sigma, so the area stays 1. The mean and the variance read off the pdf by integration are mu and sigma^2, and the width at half the peak is proportional to sigma."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-param] mu shifts the pdf, sigma stretches it and divides its \"\n      \"peak\")\neta = np.linspace(-10.0, 14.0, 240_001)\ncheck(\"adding mu shifts the pdf without changing its shape\",\n      np.allclose(pdf_gen(eta, mu, 1.0), pdf_std(eta - mu), atol=1e-12))\nfor sg in (0.5, 1.5, 3.0):\n    check(f\"sigma = {sg}: the peak is the standard peak divided by sigma\",\n          abs(pdf_gen(mu, mu, sg) - pdf_std(0.0) / sg) < 1e-12)\n    # the stretched pdf needs a correspondingly wider integration range:\n    # the fixed grid above reaches only mu + 4 sigma for sigma = 3\n    wide = np.linspace(mu - 12 * sg, mu + 12 * sg, 240_001)\n    check(f\"sigma = {sg}: the area under the pdf is 1\",\n          abs(np.trapezoid(pdf_gen(wide, mu, sg), wide) - 1.0) < 1e-9)\n\n\ndef half_max_width(sg):\n    \"\"\"Width of the pdf at half its peak, measured on the grid.\"\"\"\n    dens = pdf_gen(eta, mu, sg)\n    above = eta[dens >= dens.max() / 2.0]\n    return above[-1] - above[0]\n\n\nw1 = half_max_width(1.0)\ncheck(\"the width at half the peak is proportional to sigma\",\n      all(abs(half_max_width(sg) - sg * w1) < 1e-3 for sg in (0.5, 1.5, 3.0)))\ndens = pdf_gen(eta, mu, sigma)\ncheck(\"integrating eta against the pdf returns the mean mu\",\n      abs(np.trapezoid(eta * dens, eta) - mu) < 1e-6)\ncheck(\"integrating (eta - mu)^2 against the pdf returns sigma^2\",\n      abs(np.trapezoid((eta - mu) ** 2 * dens, eta) - sigma**2) < 1e-6)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-vec]** A Gaussian random vector x := A z + mu, with z a vector of i.i.d. standard Gaussian RVs and A a matrix square root of C (a Cholesky factor), has mean mu and covariance matrix C. Its components are independent exactly when C is diagonal: for a diagonal C the joint probabilities factorize into the marginals, for a non-diagonal one they do not."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-vec] x := A z + mu has covariance matrix C = A A^T; \"\n      \"independent components exactly for a diagonal C\")\nC = np.array([[1.0, 0.6, 0.2], [0.6, 2.0, -0.3], [0.2, -0.3, 0.8]])\nA = np.linalg.cholesky(C)\nmuvec = np.array([1.0, -2.0, 0.5])\ncheck(\"A is a matrix square root of C: A A^T = C\",\n      np.allclose(A @ A.T, C, atol=1e-12))\nzv = rng.standard_normal((10**6, 3))\nxv = zv @ A.T + muvec\ncheck(\"the sample mean of x matches mu\",\n      np.max(np.abs(xv.mean(axis=0) - muvec)) < 1e-2)\nxc = xv - xv.mean(axis=0)\nC_emp = xc.T @ xc / xv.shape[0]\ncheck(\"the sample covariance matrix matches C\",\n      np.max(np.abs(C_emp - C)) < 2e-2)\ncorr12 = C_emp[0, 1] / np.sqrt(C_emp[0, 0] * C_emp[1, 1])\ncheck(\"for this non-diagonal C the first two components are correlated, \"\n      \"hence dependent\", corr12 > 0.3)\n\n\ndef factorization_gap(cov):\n    \"\"\"Largest gap between the joint probabilities of two components and\n    the product of their separate probabilities, on a coarse grid.\"\"\"\n    w = rng.standard_normal((400_000, 2)) @ np.linalg.cholesky(cov).T\n    edges = [np.quantile(w[:, j], np.linspace(0.0, 1.0, 7)) for j in (0, 1)]\n    edges[0][0] = edges[1][0] = -np.inf\n    edges[0][-1] = edges[1][-1] = np.inf\n    i0 = np.digitize(w[:, 0], edges[0][1:-1])\n    i1 = np.digitize(w[:, 1], edges[1][1:-1])\n    joint = np.zeros((6, 6))\n    np.add.at(joint, (i0, i1), 1.0 / w.shape[0])\n    return float(np.max(np.abs(joint - np.outer(joint.sum(1), joint.sum(0)))))\n\n\ngap_diag = factorization_gap(np.diag([1.0, 2.0]))\ngap_corr = factorization_gap(np.array([[1.0, 0.8], [0.8, 2.0]]))\ncheck(\"for a diagonal C the joint probabilities factorize \"\n      f\"(largest gap {gap_diag:.4f})\", gap_diag < 5e-3)\ncheck(\"for a non-diagonal C they do not \"\n      f\"(largest gap {gap_corr:.4f}, two orders of magnitude larger)\",\n      gap_corr > 2e-2)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-gp]** A Gaussian random vector is a stochastic process indexed by {1, ..., d}: its restriction to a subset of the indices is again a Gaussian random vector, with the covariance matrix cut out of C along those indices."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-gp] the restriction to a subset of the indices is again a \"\n      \"Gaussian random vector\")\nsub = [0, 2]\nC_sub = C[np.ix_(sub, sub)]\nxs = xv[:, sub]\nxs_c = xs - xs.mean(axis=0)\ncheck(\"the sample covariance of components 1 and 3 matches the \"\n      \"corresponding part of C\",\n      np.max(np.abs(xs_c.T @ xs_c / xs.shape[0] - C_sub)) < 2e-2)\ncheck(\"the fraction of realizations of component 3 below 1.0 matches \"\n      \"Phi((1.0 - mu_3)/sqrt(C_33))\",\n      abs(np.mean(xv[:, 2] <= 1.0)\n          - float(Phi((1.0 - muvec[2]) / np.sqrt(C[2, 2])))) < 3e-3)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-conc]** A Lipschitz function of a standard Gaussian random vector concentrates: the Euclidean norm has Lipschitz constant 1, so the length of z in d dimensions has expectation at most sqrt(d), variance at most 1, and leaves a window of width t around its expectation with probability at most 2 exp(-t^2/2) \u2014 at d = 10, 100 and 1000 alike. The squared length, which is not Lipschitz, instead fluctuates ten times more at d = 1000 than at d = 10. The bound constrains the tails only: the largest of 50 standard Gaussian RVs obeys it and is still skewed, hence not Gaussian."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-conc] a Lipschitz function of a standard Gaussian random \"\n      \"vector stays near its expectation, with a bound free of d\")\npair_a = rng.standard_normal((50_000, 20))\npair_b = rng.standard_normal((50_000, 20))\ncheck(\"the Euclidean norm is Lipschitz with constant 1\",\n      np.all(np.abs(np.linalg.norm(pair_a, axis=1)\n                    - np.linalg.norm(pair_b, axis=1))\n             <= np.linalg.norm(pair_a - pair_b, axis=1) + 1e-12))\nsd_len, sd_sqlen, dims = [], [], (10, 100, 1000)\nfor dim, m_dim in zip(dims, (200_000, 100_000, 20_000)):\n    zd = rng.standard_normal((m_dim, dim))\n    length = np.linalg.norm(zd, axis=1)\n    check(f\"d = {dim}: the expectation of the length is at most sqrt(d)\",\n          length.mean() <= np.sqrt(dim))\n    check(f\"d = {dim}: its variance stays below the squared Lipschitz \"\n          \"constant 1\", np.var(length) < 1.0)\n    for t in (1.0, 2.0, 3.0):\n        check(f\"d = {dim}, t = {t}: the fraction of realizations further \"\n              \"than t from the mean is below 2 exp(-t^2/2)\",\n              np.mean(np.abs(length - length.mean()) >= t)\n              <= 2.0 * np.exp(-(t**2) / 2.0))\n    sd_len.append(float(np.std(length)))\n    sd_sqlen.append(float(np.std(np.sum(zd**2, axis=1))))\ncheck(\"the fluctuation of the length is the same at every d \"\n      f\"(standard deviations {', '.join(f'{v:.2f}' for v in sd_len)})\",\n      max(sd_len) - min(sd_len) < 0.02)\ncheck(\"the squared length is not Lipschitz and fluctuates ten times \"\n      f\"more at d = 1000 than at d = 10 (factor \"\n      f\"{sd_sqlen[2] / sd_sqlen[0]:.1f})\",\n      9.0 < sd_sqlen[2] / sd_sqlen[0] < 11.0)\nlargest = np.max(rng.standard_normal((200_000, 50)), axis=1)\ncheck(\"the largest of 50 standard Gaussian RVs obeys the same bound\",\n      all(np.mean(np.abs(largest - largest.mean()) >= t)\n          <= 2.0 * np.exp(-(t**2) / 2.0) for t in (1.0, 2.0, 3.0)))\nskewness = float(np.mean((largest - largest.mean()) ** 3) / largest.std() ** 3)\ncheck(f\"yet its distribution is skewed ({skewness:.2f}), so the bound \"\n      \"does not make it Gaussian\", skewness > 0.3)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-ent]** Among RVs with a given variance, the Gaussian one maximizes the differential entropy: by numerical integration its entropy is log(2 pi e sigma^2)/2, above that of a uniform RV and of a sum of two uniform RVs with the same variance."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-ent] among RVs with a given variance, the Gaussian one \"\n      \"maximizes the differential entropy\")\nvar_fixed = 1.0\n\n\ndef diff_entropy(dens_vals, pts):\n    \"\"\"-integral p log p, by numerical integration.\"\"\"\n    safe = np.where(dens_vals > 0.0, dens_vals, 1.0)\n    return float(-np.trapezoid(dens_vals * np.log(safe), pts))\n\n\npts = np.linspace(-12.0, 12.0, 480_001)\nh_gauss = diff_entropy(pdf_gen(pts, 0.0, np.sqrt(var_fixed)), pts)\ncheck(\"the Gaussian differential entropy is log(2 pi e sigma^2)/2\",\n      abs(h_gauss - 0.5 * np.log(2.0 * np.pi * np.e * var_fixed)) < 1e-6)\nhalf = np.sqrt(3.0 * var_fixed)                  # uniform on [-half, half]\ndens_unif = np.where(np.abs(pts) <= half, 1.0 / (2.0 * half), 0.0)\n# the grid does not fall on +-sqrt(3), where this density jumps, so the\n# numerical integration of a step function is accurate to ~1e-4, not 1e-6\ncheck(\"that uniform RV has variance 1\",\n      abs(np.trapezoid(pts**2 * dens_unif, pts) - var_fixed) < 1e-3)\nh_unif = diff_entropy(dens_unif, pts)\na = np.sqrt(1.5 * var_fixed)                     # sum of two U[-a, a]\ndens_sum = np.where(np.abs(pts) <= 2 * a,\n                    (2 * a - np.abs(pts)) / (4 * a**2), 0.0)\ncheck(\"that sum of two uniform RVs has variance 1\",\n      abs(np.trapezoid(pts**2 * dens_sum, pts) - var_fixed) < 1e-6)\nh_sum = diff_entropy(dens_sum, pts)\ncheck(\"the Gaussian entropy exceeds the uniform one\", h_gauss > h_unif)\ncheck(\"the Gaussian entropy exceeds that of the sum of two uniform RVs\",\n      h_gauss > h_sum)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-reg]** For the regression model y = w^T x + eps with eps a Gaussian RV of variance sigma^2, the negative log-likelihood of the labels of the training set is an increasing affine function of the average squared error loss. Maximizing the likelihood over w therefore returns the same weights as minimizing that loss, which is the ERM problem of linear regression."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-reg] with Gaussian noise, maximizing the likelihood minimizes \"\n      \"the average squared error loss\")\nm_train, d_feat = 400, 3\nX = rng.standard_normal((m_train, d_feat))\nw_true = np.array([1.5, -0.7, 2.0])\nnoise_sd = 0.4\ny = X @ w_true + noise_sd * rng.standard_normal(m_train)\n\n\ndef avg_sqerr(w):\n    return float(np.mean((y - X @ w) ** 2))\n\n\ndef neg_loglik(w):\n    r = y - X @ w\n    return float(m_train / 2.0 * np.log(2.0 * np.pi * noise_sd**2)\n                 + np.sum(r**2) / (2.0 * noise_sd**2))\n\n\nconst = m_train / 2.0 * np.log(2.0 * np.pi * noise_sd**2)\nslope = m_train / (2.0 * noise_sd**2)\nprobe = [w_true, np.zeros(d_feat), np.array([0.5, 0.5, 0.5]),\n         w_true + 0.3, rng.standard_normal(d_feat)]\ncheck(\"the negative log-likelihood is const + slope * (average squared \"\n      \"error loss)\",\n      all(abs(neg_loglik(w) - (const + slope * avg_sqerr(w))) < 1e-8\n          for w in probe))\nw_hat = np.linalg.lstsq(X, y, rcond=None)[0]\ncheck(\"the least-squares weights minimize the average squared error loss \"\n      \"among the probes\",\n      all(avg_sqerr(w_hat) <= avg_sqerr(w) + 1e-12 for w in probe))\ncheck(\"they also maximize the likelihood among the probes\",\n      all(neg_loglik(w_hat) <= neg_loglik(w) + 1e-8 for w in probe))\ncheck(\"perturbing them in any direction lowers the likelihood\",\n      all(neg_loglik(w_hat + 0.05 * rng.standard_normal(d_feat))\n          > neg_loglik(w_hat) for _ in range(20)))\ncheck(\"the learned weights are close to the ones used to generate the \"\n      \"labels\", np.max(np.abs(w_hat - w_true)) < 0.1)\n\n# ------------------------------------------------------------- preview\nfig, ax = plt.subplots(1, 4, figsize=(16.6, 3.4))\nax[0].bar(s16, p16 / width16, width=0.8 * width16, facecolor=\"white\",\n          edgecolor=\"black\", label=\"sum of 16 coin flips\")\nax[0].plot(np.linspace(-4, 4, 401), pdf_std(np.linspace(-4, 4, 401)),\n           \"k--\", label=\"standard Gaussian pdf\")\nax[0].set_xlabel(\"normalized sum $(2k-n)/\\\\sqrt{n}$\")\nax[0].set_ylabel(\"probability density\")\nax[0].set_title(\"[P-gen] 16 coin flips against the Gaussian pdf\",\n                fontsize=10)\nax[0].legend(frameon=False, loc=\"upper left\", fontsize=8)\nax[1].plot(n_grid, gap_grid, \"ko-\", markersize=3,\n           label=\"largest gap to the Gaussian CDF\")\nax[1].set_xlabel(\"number of coin flips $n$\")\nax[1].set_ylabel(\"largest gap between the CDFs\")\nax[1].set_title(\"[P-gen] the gap shrinks with $n$, and stays positive\",\n                fontsize=10)\nax[1].legend(frameon=False, fontsize=8)\nax[2].plot(dims, sd_len, \"ko-\", label=\"length $\\\\|z\\\\|$ (Lipschitz)\")\nax[2].plot(dims, sd_sqlen, \"ks--\", markerfacecolor=\"white\",\n           label=\"squared length (not Lipschitz)\")\nax[2].set_xscale(\"log\")\nax[2].set_yscale(\"log\")\nax[2].set_xlabel(\"number of arguments $d$\")\nax[2].set_ylabel(\"standard deviation of the value\")\nax[2].set_title(\"[P-conc] only the Lipschitz one stays flat\", fontsize=10)\nax[2].legend(frameon=False, fontsize=8)\nnames = [\"Gaussian\", \"uniform\", \"sum of two\\nuniform RVs\"]\nax[3].bar(range(3), [h_gauss, h_unif, h_sum], width=0.6, facecolor=\"white\",\n          edgecolor=\"black\", hatch=\"///\")\nfor i, h in enumerate([h_gauss, h_unif, h_sum]):\n    ax[3].text(i, h + 0.02, f\"{h:.3f}\", ha=\"center\", fontsize=8)\nax[3].set_xticks(range(3))\nax[3].set_xticklabels(names, fontsize=8)\nax[3].set_xlabel(\"RV with variance 1\")\nax[3].set_ylabel(\"differential entropy (natural log)\")\nax[3].set_ylim(0, 1.75)\nax[3].set_title(\"[P-ent] largest entropy at a fixed variance\",\n                fontsize=10)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"gaussrv.png\", dpi=110)\nprint(f\"\\n{sum(ok for _, ok in report)}/{len(report)} checks passed\")\nassert all(ok for _, ok in report)"
  }
 ]
}