Bayesian Inference

Priors, Posteriors, Shrinkage, Model Selection, and Decisions

We Have Already Used a Prior Once

// ---------------------------------------------------------------------------
// ldk — a tiny, dependency-free plotting + statistics library for these decks.
//
// Everything runs in the browser as plain JavaScript building SVG nodes, so a
// slider move costs a fraction of a millisecond instead of a round trip
// through a WebAssembly R runtime.
// ---------------------------------------------------------------------------
ldk = {
  const NS = "http://www.w3.org/2000/svg";
  let uid = 0;

  // Deck palette, matching ldkslides.css
  const C = {
    green:  "#74b83d",
    blue:   "#4072c2",
    pink:   "#ED017D",
    orange: "#d18b2a",
    grey:   "#6E6E6E",
    ink:    "#222222",
    axis:   "#9aa0a6",
    faint:  "#efefef"
  };

  function el(tag, attrs, text) {
    const e = document.createElementNS(NS, tag);
    if (attrs) for (const k in attrs) if (attrs[k] != null) e.setAttribute(k, attrs[k]);
    if (text != null)
      // SVG collapses runs of whitespace, which would eat the wide gaps that
      // separate the readout columns; non-breaking spaces survive. (xml:space
      // is a namespaced attribute and does not survive setAttribute here.)
      e.textContent = tag === "text"
        ? String(text).replace(/ {2,}/g, m => " ".repeat(m.length))
        : text;
    return e;
  }

  function alpha(hex, a) {
    const n = parseInt(hex.slice(1), 16);
    return `rgba(${(n >> 16) & 255},${(n >> 8) & 255},${n & 255},${a})`;
  }

  // ------------------------------------------------------------------ numerics

  const LGC = [676.5203681218851, -1259.1392167224028, 771.32342877765313,
               -176.61502916214059, 12.507343278686905, -0.13857109526572012,
               9.9843695780195716e-6, 1.5056327351493116e-7];

  function lgamma(z) {
    if (z < 0.5) return Math.log(Math.PI / Math.sin(Math.PI * z)) - lgamma(1 - z);
    z -= 1;
    let x = 0.99999999999980993;
    for (let i = 0; i < 8; i++) x += LGC[i] / (z + i + 1);
    const t = z + 7.5;
    return 0.5 * Math.log(2 * Math.PI) + (z + 0.5) * Math.log(t) - t + Math.log(x);
  }

  const lchoose = (n, k) => lgamma(n + 1) - lgamma(k + 1) - lgamma(n - k + 1);

  // Chebyshev erfc (Numerical Recipes); plenty of accuracy for plotting.
  const ERFC_COF = [
    -1.3026537197817094, 6.4196979235649026e-1, 1.9476473204185836e-2,
    -9.561514786808631e-3, -9.46595344482036e-4, 3.66839497852761e-4,
     4.2523324806907e-5, -2.0278578112534e-5, -1.624290004647e-6,
     1.303655835580e-6, 1.5626441722e-8, -8.5238095915e-8,
     6.529054439e-9, 5.059343495e-9, -9.91364156e-10, -2.27365122e-10,
     9.6467911e-11, 2.394038e-12, -6.886027e-12, 8.94487e-13,
     3.13092e-13, -1.12708e-13, 3.81e-16, 7.106e-15];

  function erfc(x) {
    const z = Math.abs(x), t = 2 / (2 + z), ty = 4 * t - 2;
    let d = 0, dd = 0;
    for (let j = ERFC_COF.length - 1; j > 0; j--) {
      const tmp = d;
      d = ty * d - dd + ERFC_COF[j];
      dd = tmp;
    }
    const ans = t * Math.exp(-z * z + 0.5 * (ERFC_COF[0] + ty * d) - dd);
    return x >= 0 ? ans : 2 - ans;
  }

  const dnorm = (x, m = 0, s = 1) => {
    const z = (x - m) / s;
    return Math.exp(-0.5 * z * z) / (s * Math.sqrt(2 * Math.PI));
  };
  const pnorm = (x, m = 0, s = 1) => 0.5 * erfc(-(x - m) / (s * Math.SQRT2));

  // Acklam's inverse normal CDF (relative error < 1.2e-9)
  const QN_A = [-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02,
                 1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00];
  const QN_B = [-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02,
                 6.680131188771972e+01, -1.328068155288572e+01];
  const QN_C = [-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00,
                -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00];
  const QN_D = [7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00,
                3.754408661907416e+00];

  function qnorm(p, m = 0, s = 1) {
    if (p <= 0) return -Infinity;
    if (p >= 1) return Infinity;
    const lo = 0.02425, hi = 1 - lo;
    let x, q, r;
    if (p < lo) {
      q = Math.sqrt(-2 * Math.log(p));
      x = (((((QN_C[0] * q + QN_C[1]) * q + QN_C[2]) * q + QN_C[3]) * q + QN_C[4]) * q + QN_C[5]) /
          ((((QN_D[0] * q + QN_D[1]) * q + QN_D[2]) * q + QN_D[3]) * q + 1);
    } else if (p <= hi) {
      q = p - 0.5; r = q * q;
      x = (((((QN_A[0] * r + QN_A[1]) * r + QN_A[2]) * r + QN_A[3]) * r + QN_A[4]) * r + QN_A[5]) * q /
          (((((QN_B[0] * r + QN_B[1]) * r + QN_B[2]) * r + QN_B[3]) * r + QN_B[4]) * r + 1);
    } else {
      q = Math.sqrt(-2 * Math.log(1 - p));
      x = -(((((QN_C[0] * q + QN_C[1]) * q + QN_C[2]) * q + QN_C[3]) * q + QN_C[4]) * q + QN_C[5]) /
           ((((QN_D[0] * q + QN_D[1]) * q + QN_D[2]) * q + QN_D[3]) * q + 1);
    }
    return m + s * x;
  }

  function dbinom(k, n, p) {
    if (k < 0 || k > n) return 0;
    if (p <= 0) return k === 0 ? 1 : 0;
    if (p >= 1) return k === n ? 1 : 0;
    return Math.exp(lchoose(n, k) + k * Math.log(p) + (n - k) * Math.log1p(-p));
  }
  const dpois = (k, lam) => k < 0 ? 0 : Math.exp(-lam + k * Math.log(lam) - lgamma(k + 1));
  const dexp  = (x, rate) => x < 0 ? 0 : rate * Math.exp(-rate * x);
  const pexp  = (x, rate) => x <= 0 ? 0 : 1 - Math.exp(-rate * x);
  const dunif = (x, a, b) => (x >= a && x <= b) ? 1 / (b - a) : 0;
  const punif = (x, a, b) => x <= a ? 0 : x >= b ? 1 : (x - a) / (b - a);
  const dlnorm = (x, ml = 0, sl = 1) => x <= 0 ? 0 : dnorm(Math.log(x), ml, sl) / x;
  const dchisq1 = (x) => x <= 0 ? 0 : dnorm(Math.sqrt(x)) / Math.sqrt(x);

  // ---- regularized lower incomplete gamma P(a,x): series + continued fraction
  function gammaP(a, x) {
    if (!(x >= 0) || !(a > 0)) return NaN;
    if (x === 0) return 0;
    const pref = Math.exp(-x + a * Math.log(x) - lgamma(a));
    if (x < a + 1) {                       // series
      let ap = a, sum = 1 / a, del = sum;
      for (let i = 0; i < 500; i++) {
        ap++; del *= x / ap; sum += del;
        if (Math.abs(del) < Math.abs(sum) * 1e-16) break;
      }
      return sum * pref;
    }
    const TINY = 1e-300;                   // continued fraction for Q = 1 - P
    let b = x + 1 - a, c = 1 / TINY, d = 1 / b, h = d;
    for (let i = 1; i <= 500; i++) {
      const an = -i * (i - a);
      b += 2;
      d = an * d + b; if (Math.abs(d) < TINY) d = TINY;
      c = b + an / c; if (Math.abs(c) < TINY) c = TINY;
      d = 1 / d;
      const del = d * c;
      h *= del;
      if (Math.abs(del - 1) < 1e-16) break;
    }
    return 1 - pref * h;
  }

  // ---- regularized incomplete beta I_x(a,b)
  function betacf(a, b, x) {
    const TINY = 1e-300, qab = a + b, qap = a + 1, qam = a - 1;
    let c = 1, d = 1 - qab * x / qap;
    if (Math.abs(d) < TINY) d = TINY;
    d = 1 / d;
    let h = d;
    for (let m = 1; m <= 400; m++) {
      const m2 = 2 * m;
      let aa = m * (b - m) * x / ((qam + m2) * (a + m2));
      d = 1 + aa * d; if (Math.abs(d) < TINY) d = TINY;
      c = 1 + aa / c; if (Math.abs(c) < TINY) c = TINY;
      d = 1 / d; h *= d * c;
      aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
      d = 1 + aa * d; if (Math.abs(d) < TINY) d = TINY;
      c = 1 + aa / c; if (Math.abs(c) < TINY) c = TINY;
      d = 1 / d;
      const del = d * c;
      h *= del;
      if (Math.abs(del - 1) < 3e-16) break;
    }
    return h;
  }
  function betai(a, b, x) {
    if (!(x > 0)) return 0;
    if (x >= 1) return 1;
    const bt = Math.exp(lgamma(a + b) - lgamma(a) - lgamma(b) +
                        a * Math.log(x) + b * Math.log1p(-x));
    return x < (a + 1) / (a + b + 2)
      ? bt * betacf(a, b, x) / a
      : 1 - bt * betacf(b, a, 1 - x) / b;
  }

  // ---- chi-squared and Student t
  const dchisq = (x, k) => x <= 0 ? 0 :
    Math.exp((k / 2 - 1) * Math.log(x) - x / 2 - lgamma(k / 2) - (k / 2) * Math.LN2);
  const pchisq = (x, k) => x <= 0 ? 0 : gammaP(k / 2, x / 2);

  const dt = (x, nu) => Math.exp(lgamma((nu + 1) / 2) - lgamma(nu / 2)
    - 0.5 * Math.log(nu * Math.PI) - (nu + 1) / 2 * Math.log1p(x * x / nu));
  const pt = (x, nu) => {
    const half = 0.5 * betai(nu / 2, 0.5, nu / (nu + x * x));
    return x > 0 ? 1 - half : half;
  };

  // quantiles by bisection — a handful of calls per redraw, so speed is moot
  function invert(cdf, p, lo, hi) {
    for (let i = 0; i < 90; i++) {
      const mid = (lo + hi) / 2;
      if (cdf(mid) < p) lo = mid; else hi = mid;
    }
    return (lo + hi) / 2;
  }
  const qchisq = (p, k) => invert(x => pchisq(x, k), p, 0, Math.max(60, k * 25));
  const qt = (p, nu) => invert(x => pt(x, nu), p, -1e4, 1e4);

  // Cumulative sums of a pmf enumerated over 0..n — cheaper and more accurate
  // here than a general incomplete-beta / incomplete-gamma implementation.
  function cdfFromPmf(pmf) {
    const out = new Array(pmf.length);
    let s = 0;
    for (let i = 0; i < pmf.length; i++) { s += pmf[i]; out[i] = Math.min(s, 1); }
    return out;
  }

  // ----------------------------------------------------------------- sequences

  function seq(a, b, n) {
    const out = new Array(n);
    for (let i = 0; i < n; i++) out[i] = a + (b - a) * i / (n - 1);
    return out;
  }
  const cumsum = (v) => { let s = 0; return v.map(x => (s += x)); };
  const sum = (v) => v.reduce((a, b) => a + b, 0);
  const mean = (v) => sum(v) / v.length;

  function sd(v) {
    const m = mean(v);
    return Math.sqrt(v.reduce((a, x) => a + (x - m) * (x - m), 0) / (v.length - 1));
  }

  function quantileSorted(s, p) {
    const h = (s.length - 1) * p, i = Math.floor(h);
    const j = Math.min(i + 1, s.length - 1);
    return s[i] + (h - i) * (s[j] - s[i]);
  }

  function hist(v, lo, hi, nb) {
    const cnt = new Float64Array(nb), w = (hi - lo) / nb;
    for (let i = 0; i < v.length; i++) {
      const k = Math.floor((v[i] - lo) / w);
      if (k >= 0 && k < nb) cnt[k]++;
      else if (v[i] === hi) cnt[nb - 1]++;
    }
    const dens = new Array(nb);
    for (let i = 0; i < nb; i++) dens[i] = cnt[i] / (v.length * w);
    return { lo, hi, w, nb, dens, count: cnt };
  }

  // --------------------------------------------------------------------- rngs

  // mulberry32: small, fast, seedable — the demos stay reproducible.
  function rng(seed) {
    let a = seed >>> 0;
    const u = () => {
      a |= 0; a = (a + 0x6D2B79F5) | 0;
      let t = Math.imul(a ^ (a >>> 15), 1 | a);
      t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t;
      return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
    };
    let spare = null;
    u.normal = () => {                       // Marsaglia polar
      if (spare !== null) { const v = spare; spare = null; return v; }
      let x, y, s;
      do { x = 2 * u() - 1; y = 2 * u() - 1; s = x * x + y * y; } while (s >= 1 || s === 0);
      const f = Math.sqrt(-2 * Math.log(s) / s);
      spare = y * f;
      return x * f;
    };
    u.exp = () => -Math.log(1 - u());
    u.cauchy = () => Math.tan(Math.PI * (u() - 0.5));
    return u;
  }

  // Weighted sampling without replacement (Efraimidis–Spirakis exponential race)
  function sampleWeighted(n, w, u) {
    const N = w.length;
    const keys = new Array(N);
    for (let i = 0; i < N; i++) keys[i] = w[i] > 0 ? -Math.log(1 - u()) / w[i] : Infinity;
    const idx = Array.from({ length: N }, (_, i) => i);
    idx.sort((a, b) => keys[a] - keys[b]);
    return idx.slice(0, n);
  }

  // Sampling *with* replacement — what the bootstrap needs, and what
  // sampleWeighted deliberately is not.
  function sampleReplace(n, v, u) {
    const out = new Array(n);
    for (let i = 0; i < n; i++) out[i] = v[Math.floor(u() * v.length)];
    return out;
  }

  // Beta density — the conjugate posterior for a proportion.
  const dbeta = (x, a, b) => (x <= 0 || x >= 1) ? 0 :
    Math.exp(lgamma(a + b) - lgamma(a) - lgamma(b) +
             (a - 1) * Math.log(x) + (b - 1) * Math.log(1 - x));

  // Beta quantile, by inverting the regularized incomplete beta above. Used for
  // the Clopper-Pearson interval, whose endpoints are exactly beta quantiles.
  const qbeta = (p, a, b) => invert(x => betai(a, b, x), p, 0, 1);

  // Gini coefficient of a sorted, non-negative sample: twice the area between
  // the Lorenz curve and the diagonal. No closed-form standard error exists,
  // which is exactly why the bootstrap earns its place.
  function gini(sorted) {
    const n = sorted.length;
    let num = 0, tot = 0;
    for (let i = 0; i < n; i++) { num += (i + 1) * sorted[i]; tot += sorted[i]; }
    return tot > 0 ? (2 * num) / (n * tot) - (n + 1) / n : 0;
  }

  // ------------------------------------------------------- regression + series

  // Simple OLS of y on x, with both the classical and the heteroskedasticity-
  // robust (HC1) standard error. Deck 08 leans on the pair: the estimate and its
  // standard error fail independently, and the standard error fails first.
  function ols(x, y) {
    const n = x.length;
    const mx = mean(x), my = mean(y);
    let sxx = 0, sxy = 0;
    for (let i = 0; i < n; i++) { const dx = x[i] - mx; sxx += dx * dx; sxy += dx * (y[i] - my); }
    const b = sxx > 0 ? sxy / sxx : 0;
    const a = my - b * mx;
    let rss = 0, meat = 0, tss = 0;
    for (let i = 0; i < n; i++) {
      const e = y[i] - a - b * x[i], dx = x[i] - mx;
      rss += e * e; meat += dx * dx * e * e; tss += (y[i] - my) * (y[i] - my);
    }
    const se = sxx > 0 && n > 2 ? Math.sqrt(rss / (n - 2) / sxx) : NaN;
    const seR = sxx > 0 && n > 2
      ? Math.sqrt(meat * n / (n - 2)) / sxx                    // HC1
      : NaN;
    return { a, b, se, seR, t: b / se, tR: b / seR,
             r2: tss > 0 ? 1 - rss / tss : 0, n };
  }

  const corr = (x, y) => {
    const mx = mean(x), my = mean(y);
    let sxy = 0, sxx = 0, syy = 0;
    for (let i = 0; i < x.length; i++) {
      const dx = x[i] - mx, dy = y[i] - my;
      sxy += dx * dy; sxx += dx * dx; syy += dy * dy;
    }
    return sxx > 0 && syy > 0 ? sxy / Math.sqrt(sxx * syy) : 0;
  };

  // Stationary AR(1), started from its stationary distribution so the first
  // observations are not systematically closer to zero than the rest.
  function ar1(n, rho, u, sd = 1) {
    const out = new Float64Array(n);
    let v = u.normal() * sd / Math.sqrt(Math.max(1 - rho * rho, 1e-9));
    for (let i = 0; i < n; i++) { v = rho * v + sd * u.normal(); out[i] = v; }
    return out;
  }

  const randomWalk = (n, u, sd = 1) => {
    const out = new Float64Array(n);
    let v = 0;
    for (let i = 0; i < n; i++) { v += sd * u.normal(); out[i] = v; }
    return out;
  };

  // G clusters of m units with intra-cluster correlation rho: a shared group
  // draw plus an idiosyncratic one, weighted so the total variance stays 1.
  function clustered(G, m, rho, u) {
    const out = new Float64Array(G * m);
    const sg = Math.sqrt(Math.max(rho, 0)), si = Math.sqrt(Math.max(1 - rho, 0));
    for (let g = 0; g < G; g++) {
      const shared = sg * u.normal();
      for (let i = 0; i < m; i++) out[g * m + i] = shared + si * u.normal();
    }
    return out;
  }

  // ---------------------------------------------------------------- formatting

  const fmt = (v, d = 2) => (Number.isFinite(v) ? v.toFixed(d) : "—");
  const signed = (v, d = 2) => (Number.isFinite(v) ? (v >= 0 ? "+" : "") + v.toFixed(d) : "—");

  function smart(v) {
    if (!Number.isFinite(v)) return "—";
    const a = Math.abs(v);
    if (a >= 1e5 || (a > 0 && a < 1e-3)) return v.toExponential(1);
    return v.toFixed(a >= 100 ? 0 : a >= 10 ? 1 : 2);
  }

  function nice(lo, hi, n) {
    const span = (hi - lo) || Math.abs(hi) || 1;
    const raw = span / n;
    const mag = Math.pow(10, Math.floor(Math.log10(raw)));
    const norm = raw / mag;
    const step = mag * (norm < 1.5 ? 1 : norm < 3 ? 2 : norm < 7 ? 5 : 10);
    const out = [];
    for (let v = Math.ceil(lo / step - 1e-9) * step; v <= hi + 1e-9 * span; v += step)
      out.push(Math.abs(v) < step * 1e-9 ? 0 : v);
    return { ticks: out, step };
  }

  function logTicks(lo, hi) {
    const decades = Math.log10(hi / lo);
    const mults = decades > 2.2 ? [1] : decades > 1.1 ? [1, 3] : [1, 2, 5];
    const out = [];
    for (let e = Math.floor(Math.log10(lo)); e <= Math.ceil(Math.log10(hi)); e++)
      for (const k of mults) {
        const v = k * Math.pow(10, e);
        if (v >= lo * 0.9999 && v <= hi * 1.0001) out.push(v);
      }
    return out;
  }

  const decimalsFor = (step) =>
    Math.max(0, Math.min(6, -Math.floor(Math.log10(Math.abs(step)) + 1e-9)));

  // -------------------------------------------------------------------- charts

  function chart(W, H, extra) {
    const svg = el("svg", Object.assign({
      viewBox: `0 0 ${W} ${H}`,
      style: `width:100%;height:auto;display:block;font-family:inherit;color:${C.ink}`
    }, extra || {}));
    svg.__W = W;
    svg.__H = H;
    return svg;
  }

  function panel(svg, o) {
    const L = o.left ?? 0, T = o.top ?? 0;
    const W = o.w ?? svg.__W, H = o.h ?? svg.__H;
    const m = Object.assign({ l: 74, r: 22, t: 52, b: 54 }, o.margin || {});
    const fs = o.fontSize ?? 15;
    const iw = W - m.l - m.r, ih = H - m.t - m.b;
    const px0 = L + m.l, py0 = T + m.t;

    let [xa, xb] = o.xdom, [ya, yb] = o.ydom;
    if (!(xb > xa)) xb = xa + 1;
    if (!(yb > ya)) yb = ya + 1;

    const lg = o.xlog ? Math.log : (v => v);
    const xA = lg(xa), xB = lg(xb);
    const x = v => px0 + (lg(v) - xA) / (xB - xA) * iw;
    const y = v => py0 + ih - (v - ya) / (yb - ya) * ih;

    const gAxis = el("g");
    const id = `ldkclip${++uid}`;
    const defs = el("defs");
    const cp = el("clipPath", { id });
    cp.append(el("rect", { x: px0 - 3, y: py0 - 8, width: iw + 6, height: ih + 11 }));
    defs.append(cp);
    const gc = el("g", { "clip-path": `url(#${id})` });
    const gt = el("g");
    svg.append(defs, gAxis, gc, gt);

    // ticks
    let xt, xfmt, yt, yfmt;
    if (o.xticks) {
      xt = o.xticks;
      xfmt = o.xfmt || (v => String(v));
    } else if (o.xlog) {
      xt = logTicks(xa, xb);
      xfmt = o.xfmt || (v => v >= 1 ? String(Math.round(v)) : String(v));
    } else {
      const n = nice(xa, xb, o.nx ?? 6);
      xt = n.ticks;
      xfmt = o.xfmt || (v => v.toFixed(decimalsFor(n.step)));
    }
    if (o.yticks) {
      yt = o.yticks;
      yfmt = o.yfmt || (v => String(v));
    } else {
      const n = nice(ya, yb, o.ny ?? 5);
      yt = n.ticks;
      yfmt = o.yfmt || (v => v.toFixed(decimalsFor(n.step)));
    }

    if (o.axes === false) { xt = []; yt = []; }   // bare drawing surface

    if (o.axes !== false) {
    gAxis.append(el("line", { x1: px0, y1: py0 + ih, x2: px0 + iw, y2: py0 + ih,
                              stroke: C.axis, "stroke-width": 1.3 }));
    gAxis.append(el("line", { x1: px0, y1: py0, x2: px0, y2: py0 + ih,
                              stroke: C.axis, "stroke-width": 1.3 }));
    }
    for (const t of xt) {
      const px = x(t);
      if (px < px0 - 0.5 || px > px0 + iw + 0.5) continue;
      gAxis.append(el("line", { x1: px, y1: py0 + ih, x2: px, y2: py0 + ih + 5,
                                stroke: C.axis, "stroke-width": 1.3 }));
      gAxis.append(el("text", { x: px, y: py0 + ih + 6 + fs, "text-anchor": "middle",
                                "font-size": fs, fill: C.grey }, xfmt(t)));
    }
    for (const t of yt) {
      const py = y(t);
      if (py < py0 - 0.5 || py > py0 + ih + 0.5) continue;
      gAxis.append(el("line", { x1: px0 - 5, y1: py, x2: px0, y2: py,
                                stroke: C.axis, "stroke-width": 1.3 }));
      gAxis.append(el("text", { x: px0 - 10, y: py + fs * 0.35, "text-anchor": "end",
                                "font-size": fs, fill: C.grey }, yfmt(t)));
    }
    if (o.xlab)
      gAxis.append(el("text", { x: px0 + iw / 2, y: T + H - 8, "text-anchor": "middle",
                                "font-size": fs + 2, fill: C.ink }, o.xlab));
    if (o.ylab) {
      const ty = py0 + ih / 2, tx = L + 20;
      gAxis.append(el("text", { x: tx, y: ty, "text-anchor": "middle", "font-size": fs + 2,
                                fill: C.ink, transform: `rotate(-90 ${tx} ${ty})` }, o.ylab));
    }

    const P = {
      x, y, svg, gAxis, gc, gt, fs,
      xdom: [xa, xb], ydom: [ya, yb],
      area: { x0: px0, y0: py0, w: iw, h: ih },

      add(node, top) { (top ? gt : gc).append(node); return node; },

      path(pts, a) {
        let d = "";
        for (let i = 0; i < pts.length; i++)
          d += (i ? "L" : "M") + x(pts[i][0]).toFixed(2) + "," + y(pts[i][1]).toFixed(2);
        return P.add(el("path", Object.assign(
          { d, fill: "none", "stroke-linejoin": "round", "stroke-linecap": "round" }, a)));
      },

      fill(pts, base, a) {
        let d = "";
        for (let i = 0; i < pts.length; i++)
          d += (i ? "L" : "M") + x(pts[i][0]).toFixed(2) + "," + y(pts[i][1]).toFixed(2);
        d += "L" + x(pts[pts.length - 1][0]).toFixed(2) + "," + y(base).toFixed(2);
        d += "L" + x(pts[0][0]).toFixed(2) + "," + y(base).toFixed(2) + "Z";
        return P.add(el("path", Object.assign({ d, stroke: "none" }, a)));
      },

      seg(x1, y1, x2, y2, a) {
        return P.add(el("line", Object.assign(
          { x1: x(x1), y1: y(y1), x2: x(x2), y2: y(y2) }, a)));
      },
      vline(v, a) { return P.seg(v, ya, v, yb, a); },
      hline(v, a) { return P.seg(xa, v, xb, v, a); },

      bar(x1, x2, y1, y2, a) {
        const X1 = x(x1), X2 = x(x2), Y1 = y(y1), Y2 = y(y2);
        return P.add(el("rect", Object.assign({
          x: Math.min(X1, X2), y: Math.min(Y1, Y2),
          width: Math.max(Math.abs(X2 - X1), 0.5),
          height: Math.max(Math.abs(Y2 - Y1), 0)
        }, a)));
      },

      dot(vx, vy, r, a) {
        return P.add(el("circle", Object.assign({ cx: x(vx), cy: y(vy), r }, a)));
      },

      // Dirac point mass at vx drawn as an upward arrow; its height is the mass
      spike(vx, mass, a = {}) {
        const X = x(vx), Y0 = y(0), Y1 = y(mass), hw = a.head ?? 8, hl = hw * 1.7;
        const col = a.color ?? C.ink;
        const g = el("g");
        g.append(el("line", { x1: X, y1: Y0, x2: X, y2: Y1 + hl * 0.6, stroke: col,
                              "stroke-width": a.lw ?? 4, "stroke-dasharray": a.dash ?? null }));
        g.append(el("path", { d: `M${X - hw},${Y1 + hl} L${X},${Y1} L${X + hw},${Y1 + hl}Z`,
                              fill: col }));
        return P.add(g);
      },

      // mtext-style annotation above the plot region
      header(s, a = {}) {
        const adj = a.adj ?? 0;
        const py = py0 - 14 - (a.line ?? 0) * (fs + 7);
        return P.add(el("text", {
          x: adj === 1 ? px0 + iw : adj === 0.5 ? px0 + iw / 2 : px0,
          y: py,
          "text-anchor": adj === 1 ? "end" : adj === 0.5 ? "middle" : "start",
          "font-size": a.size ?? fs + 2,
          "font-weight": a.weight ?? "bold",
          fill: a.fill ?? C.ink
        }, s), true);
      },

      legend(items, opt = {}) {
        const size = opt.size ?? fs;
        const pad = opt.pad ?? 8;
        const rowH = size + 8, swW = 30;
        const g = el("g");
        let maxw = 0;
        items.forEach((it, i) => {
          const ty = i * rowH + size;
          if (it.swatch)
            g.append(el("rect", { x: 0, y: ty - size * 0.82, width: swW,
                                  height: size * 0.85, fill: it.swatch, stroke: "none" }));
          else
            g.append(el("line", { x1: 0, y1: ty - size * 0.35, x2: swW, y2: ty - size * 0.35,
                                  stroke: it.color, "stroke-width": it.lw ?? 3.5,
                                  "stroke-dasharray": it.dash ?? null,
                                  "stroke-linecap": "round" }));
          g.append(el("text", { x: swW + 9, y: ty, "font-size": size,
                                fill: opt.fill ?? C.ink }, it.label));
          maxw = Math.max(maxw, swW + 9 + it.label.length * size * 0.54);
        });
        // opt-in backing panel, for legends that must sit over the data
        if (opt.bg) {
          const q = 8;
          g.insertBefore(el("rect", {
            x: -q, y: -q + 3, width: maxw + 2 * q,
            height: items.length * rowH + 2 * q - 6,
            fill: "white", "fill-opacity": opt.bg === true ? 0.85 : opt.bg, rx: 5
          }), g.firstChild);
        }
        const corner = opt.corner ?? "topright";
        const gx = corner.includes("right") ? px0 + iw - maxw - pad : px0 + pad;
        const gy = corner.includes("top") ? py0 + pad
                                          : py0 + ih - items.length * rowH - pad;
        g.setAttribute("transform", `translate(${gx},${gy})`);
        gt.append(g);
        return g;
      }
    };
    return P;
  }

  // Convenience: one chart with exactly one panel filling it.
  function single(W, H, o) {
    const svg = chart(W, H);
    return { svg, p: panel(svg, o) };
  }

  return {
    C, el, alpha, chart, panel, single,
    lgamma, lchoose, erfc, dnorm, pnorm, qnorm, dbinom, dpois, dexp, pexp,
    dunif, punif, dlnorm, dchisq1, cdfFromPmf,
    gammaP, betai, dbeta, qbeta, dchisq, pchisq, qchisq, dt, pt, qt, invert,
    gini,
    seq, cumsum, sum, mean, sd, quantileSorted, hist, rng, sampleWeighted,
    sampleReplace, ols, corr, ar1, randomWalk, clustered,
    fmt, signed, smart, nice, logTicks
  };
}

Deck 04 shrank \(\bar X\) toward zero, \(\hat\mu_2 = c\bar X\), and showed that a little bias can lower the MSE — the idea behind ridge regression. That estimator is exactly maximum a posteriori estimation under a Gaussian prior centred at zero, and the deck never said so.

\[\hat\theta_{\text{MAP}} \;=\; \arg\max_\theta \; \underbrace{p(y \mid \theta)}_{\text{likelihood}}\;\underbrace{p(\theta)}_{\text{prior}}\]

What that slide did not say

The MAP estimate is one number read off a whole distribution. That distribution — the posterior — carries the uncertainty as well, and it answers questions a point estimate cannot.

Regularization was the posterior’s mode. This deck keeps the rest of it.

Bayes’ Theorem for Parameters

Definition — Posterior Distribution

For data \(y\) and parameter \(\theta\), \[p(\theta \mid y) \;=\; \frac{p(y \mid \theta)\,p(\theta)}{p(y)}, \qquad p(y) = \int p(y\mid\theta)\,p(\theta)\,d\theta\] The denominator does not involve \(\theta\), so in practice \[p(\theta \mid y) \;\propto\; p(y \mid \theta)\,p(\theta).\]

In words

Prior belief about \(\theta\), reweighted by how well each \(\theta\) predicted what you actually saw.

The shift

\(\theta\) is now a random variable — not because it moves, but because we are uncertain about it. Probability describes belief, not frequency.

Conjugate Updating

The Beta–Bernoulli Model

Estimating a proportion — a take-up rate, a job-finding rate — from yes/no observations.

Data

\[x_i \mid \theta \;\sim\; \text{Bernoulli}(\theta), \quad i = 1,\dots,n\] with \(k = \sum_i x_i\) successes

Prior

\[\theta \;\sim\; \text{Beta}(\alpha, \beta)\] \(\alpha, \beta > 0\) set the prior’s location and strength

Multiply the two, keeping only the parts that involve \(\theta\):

\[\underbrace{\theta^{\alpha-1}(1-\theta)^{\beta-1}}_{\text{prior}} \;\times\; \underbrace{\theta^{k}(1-\theta)^{n-k}}_{\text{likelihood}} \;=\; \theta^{\,\textcolor[RGB]{237,1,125}{\alpha+k}\,-1}\,(1-\theta)^{\,\textcolor[RGB]{237,1,125}{\beta+n-k}\,-1}\]

Posterior — the same family, with updated parameters

\[\theta \mid x \;\sim\; \text{Beta}\big(\textcolor[RGB]{237,1,125}{\alpha + k},\; \textcolor[RGB]{237,1,125}{\beta + n - k}\big)\] Successes are added to \(\alpha\), failures to \(\beta\). Nothing else about the data matters.

Beta–Bernoulli, Step by Step

Same data for both rows — 7 successes in 10 — and two different priors.

{
  const C = ldk.C;
  const n = 10, k = 7;
  const rows = [{ a: 2, b: 2, lab: "weak prior" }, { a: 20, b: 20, lab: "strong prior" }];
  const g = ldk.seq(0.002, 0.998, 400);
  const svg = ldk.chart(1180, 360);
  const W = 1180 / 3, H = 180;

  rows.forEach((r, i) => {
    const top = i * H, last = i === rows.length - 1;
    const prior = g.map(t => ldk.dbeta(t, r.a, r.b));
    const lik = g.map(t => ldk.dbeta(t, k + 1, n - k + 1));       // normalised likelihood
    const post = g.map(t => ldk.dbeta(t, r.a + k, r.b + n - k));
    const ymax = 1.12 * Math.max(...prior, ...lik, ...post);
    const mk = col => ldk.panel(svg, {
      left: col * W, top, w: W, h: H, xdom: [0, 1], ydom: [0, ymax],
      xlab: last ? "θ" : null, yticks: [], nx: 5,
      margin: { l: 30, r: 26, t: 40, b: last ? 48 : 26 }
    });
    const p1 = mk(0), p2 = mk(1), p3 = mk(2);

    p1.fill(g.map((t, j) => [t, prior[j]]), 0, { fill: ldk.alpha(C.blue, 0.18) });
    p1.path(g.map((t, j) => [t, prior[j]]), { stroke: C.blue, "stroke-width": 3 });
    p1.header(`${r.lab}:  Beta(${r.a}, ${r.b})`, { size: 16, fill: C.blue });

    p2.fill(g.map((t, j) => [t, lik[j]]), 0, { fill: ldk.alpha(C.orange, 0.18) });
    p2.path(g.map((t, j) => [t, lik[j]]), { stroke: C.orange, "stroke-width": 3 });
    p2.header(`likelihood:  ${k} of ${n}`, { size: 16, fill: C.orange });

    p3.path(g.map((t, j) => [t, prior[j]]),
            { stroke: ldk.alpha(C.blue, 0.55), "stroke-width": 2, "stroke-dasharray": "6,5" });
    p3.path(g.map((t, j) => [t, lik[j]]),
            { stroke: ldk.alpha(C.orange, 0.55), "stroke-width": 2, "stroke-dasharray": "6,5" });
    p3.fill(g.map((t, j) => [t, post[j]]), 0, { fill: ldk.alpha(C.pink, 0.2) });
    p3.path(g.map((t, j) => [t, post[j]]), { stroke: C.pink, "stroke-width": 3.4 });
    p3.header(`posterior:  Beta(${r.a + k}, ${r.b + n - k})   mean ${ldk.fmt((r.a + k) / (r.a + r.b + n), 3)}`,
              { size: 16, fill: C.pink });

    for (const [col, sym] of [[1, "×"], [2, "="]])
      svg.append(ldk.el("text", { x: col * W + 4, y: top + H / 2 + 4, "font-size": 30,
                                  "font-weight": "bold", fill: C.grey, "text-anchor": "middle" }, sym));
  });
  return svg;
}

\(\text{Beta}(2, 2) \;\to\; \text{Beta}(2 + \textcolor[RGB]{237,1,125}{7},\; 2 + \textcolor[RGB]{237,1,125}{3}) = \text{Beta}(\textcolor[RGB]{237,1,125}{9}, \textcolor[RGB]{237,1,125}{5})\)

A prior worth 4 observations against 10 real ones: the data move the mean from 0.50 to 0.64.

\(\text{Beta}(20, 20) \;\to\; \text{Beta}(20 + \textcolor[RGB]{237,1,125}{7},\; 20 + \textcolor[RGB]{237,1,125}{3}) = \text{Beta}(\textcolor[RGB]{237,1,125}{27}, \textcolor[RGB]{237,1,125}{23})\)

A prior worth 40 observations: the same data move the mean only to 0.54.

Conjugate Updating, Live

viewof cj_prior = Inputs.select(
  new Map([["Uniform — Beta(1,1)", "flat"],
           ["Weakly sceptical — Beta(2,8)", "scept"],
           ["Optimistic — Beta(8,2)", "opt"],
           ["Strong — Beta(40,40)", "strong"]]),
  {label: "Prior", value: "flat"})
viewof cj_n = Inputs.range([0, 400], {step: 1, value: 20, label: "Observations n"})
viewof cj_p = Inputs.range([0, 1], {step: 0.01, value: 0.3, label: "Observed share k/n"})
viewof cj_all = Inputs.toggle({label: "Overlay all four priors", value: false})

At \(n = 0\) the posterior is the prior. As \(n\) grows the likelihood takes over and the four posteriors converge — the data outvote the prior.

Turn the overlay on and slide \(n\) from 10 to 400 to watch it happen.

{
  const C = ldk.C;
  const n = cj_n, k = Math.round(cj_p * n);
  const priors = { flat: [1, 1], scept: [2, 8], opt: [8, 2], strong: [40, 40] };
  const names  = { flat: "Beta(1,1)", scept: "Beta(2,8)",
                   opt: "Beta(8,2)", strong: "Beta(40,40)" };
  const cols   = { flat: C.blue, scept: C.pink, opt: C.green, strong: C.orange };

  const g = ldk.seq(0.001, 0.999, 500);
  const shown = cj_all ? Object.keys(priors) : [cj_prior];

  // scaled likelihood, for display only
  const lik = g.map(t => Math.exp(k * Math.log(t) + (n - k) * Math.log(1 - t)));
  const lmax = Math.max(...lik);

  let ymax = 0;
  const post = {};
  for (const key of shown) {
    const [a, b] = priors[key];
    post[key] = g.map(t => ldk.dbeta(t, a + k, b + n - k));
    ymax = Math.max(ymax, ...post[key]);
    if (!cj_all) ymax = Math.max(ymax, ...g.map(t => ldk.dbeta(t, a, b)));
  }
  ymax = Math.min(ymax, 40);

  const { svg, p } = ldk.single(880, 520, {
    xdom: [0, 1], ydom: [0, ymax * 1.12],
    xlab: "θ", ylab: "density", margin: { t: 74 }, nx: 6
  });

  if (n > 0)
    p.fill(g.map((t, i) => [t, lik[i] / lmax * ymax * 0.9]), 0,
           { fill: ldk.alpha(C.grey, 0.16) });

  const items = [];
  for (const key of shown) {
    const [a, b] = priors[key];
    if (!cj_all)
      p.path(g.map(t => [t, ldk.dbeta(t, a, b)]),
             { stroke: cols[key], "stroke-width": 2.4, "stroke-dasharray": "7,5" });
    p.path(g.map((t, i) => [t, post[key][i]]),
           { stroke: cols[key], "stroke-width": 3.4 });
    const mean = (a + k) / (a + b + n);
    items.push({ label: `${names[key]} → mean ${ldk.fmt(mean, 3)}`,
                 color: cols[key], lw: 3.4 });
  }
  if (n > 0) items.push({ label: "scaled likelihood", swatch: ldk.alpha(C.grey, 0.16) });
  p.legend(items, { size: 14, bg: true });

  const [a0, b0] = priors[cj_prior];
  p.header(`prior ${names[cj_prior]}  +  data (${k} of ${n})` +
           `  →  posterior Beta(${a0 + k}, ${b0 + n - k})`,
           { line: 1, size: 18, fill: C.pink });
  p.header(n === 0 ? "no data yet — the posterior is the prior"
                   : `prior is worth ${a0 + b0} observations; you have ${n}`,
           { size: 15, fill: C.grey });
  return svg;
}

What the Update Does

With \(n = 0\) the posterior is exactly the prior — Bayes’ rule degrades gracefully to “we knew nothing new”
Each observation shifts \(\alpha\) or \(\beta\) by one. The prior is data you are pretending to have already seen: \(\text{Beta}(2,2)\) is worth four observations, \(\text{Beta}(50,50)\) a hundred
As \(n\) grows the posteriors become indistinguishable: the likelihood dominates any fixed prior

Bernstein–von Mises, informally

Under regularity conditions the posterior concentrates on the true \(\theta\) and becomes approximately normal with variance \(1/(nI(\theta))\) — the same Cramér–Rao quantity deck 04 derived. Asymptotically, Bayesian and frequentist answers agree.

Where it fails is where it matters

Small samples, many parameters, and priors that put zero mass on the truth. Convergence is no comfort when \(n\) is 40.

The Normal–Normal Model

Now a mean rather than a proportion, with the noise level \(\sigma\) known.

Data

\[x_i \mid \theta \;\sim\; N(\theta, \sigma^2), \quad i = 1,\dots,n\] summarised by \(\bar x\), with precision \(n/\sigma^2\)

Prior

\[\theta \;\sim\; N(\mu_0, \tau_0^2)\] with precision \(1/\tau_0^2\)

Both are exponentials of a quadratic in \(\theta\), so their product is too — a normal again:

\[\exp\!\Big(\!-\frac{(\theta-\mu_0)^2}{2\tau_0^2}\Big) \times \exp\!\Big(\!-\frac{n(\theta-\bar x)^2}{2\sigma^2}\Big) \;\propto\; \exp\!\Big(\!-\frac{(\theta-\textcolor[RGB]{237,1,125}{\mu_1})^2}{2\textcolor[RGB]{237,1,125}{\tau_1^2}}\Big)\]

Posterior — precisions add, means are weighted by precision

\[\theta \mid x \sim N(\textcolor[RGB]{237,1,125}{\mu_1}, \textcolor[RGB]{237,1,125}{\tau_1^2}), \qquad \textcolor[RGB]{237,1,125}{\tau_1^2} = \big(\tau_0^{-2} + n\,\sigma^{-2}\big)^{-1}, \qquad \textcolor[RGB]{237,1,125}{\mu_1} = \tau_1^2\,\big(\tau_0^{-2}\,\mu_0 + n\,\sigma^{-2}\,\bar x\big)\]

Normal–Normal, Step by Step

A training programme: \(n = 10\) participants, average earnings gain \(\bar x = 1.6\) (hundred €), \(\sigma = 1\). Two sceptical priors centred at zero effect.

{
  const C = ldk.C;
  const n = 10, xbar = 1.6, sigma = 1, mu0 = 0;
  const rows = [{ tau: 2, lab: "weak prior" }, { tau: 0.5, lab: "strong prior" }];
  const g = ldk.seq(-2.5, 3.5, 500);
  const svg = ldk.chart(1180, 340);
  const W = 1180 / 3, H = 170;
  const se = sigma / Math.sqrt(n);

  rows.forEach((r, i) => {
    const top = i * H, last = i === rows.length - 1;
    const prec = 1 / (r.tau * r.tau) + n / (sigma * sigma);
    const m1 = (mu0 / (r.tau * r.tau) + n * xbar / (sigma * sigma)) / prec, s1 = Math.sqrt(1 / prec);
    const prior = g.map(t => ldk.dnorm(t, mu0, r.tau));
    const lik = g.map(t => ldk.dnorm(t, xbar, se));
    const post = g.map(t => ldk.dnorm(t, m1, s1));
    const ymax = 1.12 * Math.max(...prior, ...lik, ...post);
    const mk = col => ldk.panel(svg, {
      left: col * W, top, w: W, h: H, xdom: [-2.5, 3.5], ydom: [0, ymax],
      xlab: last ? "θ" : null, yticks: [], nx: 6,
      margin: { l: 30, r: 26, t: 40, b: last ? 48 : 26 }
    });
    const p1 = mk(0), p2 = mk(1), p3 = mk(2);

    p1.fill(g.map((t, j) => [t, prior[j]]), 0, { fill: ldk.alpha(C.blue, 0.18) });
    p1.path(g.map((t, j) => [t, prior[j]]), { stroke: C.blue, "stroke-width": 3 });
    p1.header(`${r.lab}:  N(0, ${r.tau}²)`, { size: 16, fill: C.blue });

    p2.fill(g.map((t, j) => [t, lik[j]]), 0, { fill: ldk.alpha(C.orange, 0.18) });
    p2.path(g.map((t, j) => [t, lik[j]]), { stroke: C.orange, "stroke-width": 3 });
    p2.header(`likelihood:  x̄ = ${xbar}, se ${ldk.fmt(se, 2)}`, { size: 16, fill: C.orange });

    p3.path(g.map((t, j) => [t, prior[j]]),
            { stroke: ldk.alpha(C.blue, 0.55), "stroke-width": 2, "stroke-dasharray": "6,5" });
    p3.path(g.map((t, j) => [t, lik[j]]),
            { stroke: ldk.alpha(C.orange, 0.55), "stroke-width": 2, "stroke-dasharray": "6,5" });
    p3.fill(g.map((t, j) => [t, post[j]]), 0, { fill: ldk.alpha(C.pink, 0.2) });
    p3.path(g.map((t, j) => [t, post[j]]), { stroke: C.pink, "stroke-width": 3.4 });
    p3.header(`posterior:  N(${ldk.fmt(m1, 2)}, ${ldk.fmt(s1, 2)}²)`, { size: 16, fill: C.pink });

    for (const [col, sym] of [[1, "×"], [2, "="]])
      svg.append(ldk.el("text", { x: col * W + 4, y: top + H / 2 + 4, "font-size": 30,
                                  "font-weight": "bold", fill: C.grey, "text-anchor": "middle" }, sym));
  });
  return svg;
}

Precision: \(\tfrac{1}{4} + 10 = \textcolor[RGB]{237,1,125}{10.25}\) \(\;\Rightarrow\; \tau_1 = \textcolor[RGB]{237,1,125}{0.31}\)

Mean: \(\dfrac{0.25 \cdot 0 + 10 \cdot 1.6}{10.25} = \textcolor[RGB]{237,1,125}{1.56}\) — almost the data

Precision: \(4 + 10 = \textcolor[RGB]{237,1,125}{14}\) \(\;\Rightarrow\; \tau_1 = \textcolor[RGB]{237,1,125}{0.27}\)

Mean: \(\dfrac{4 \cdot 0 + 10 \cdot 1.6}{14} = \textcolor[RGB]{237,1,125}{1.14}\) — pulled toward zero

Beyond Conjugacy

When the Prior Is Not Conjugate

Conjugacy is luck: prior \(\times\) likelihood lands back in the prior’s family. Most pairs do not.

Example. The same yes/no data, but the prior is set on the log-odds \(\eta = \log\frac{\theta}{1-\theta}\) — natural when \(\theta\) comes out of a logit model:

\[\eta \sim N(m, s^2)\]

\[p(\eta \mid x) \;\propto\; \frac{e^{k\eta}}{(1+e^{\eta})^{n}}\; \exp\!\Big(\!-\frac{(\eta-m)^2}{2s^2}\Big)\]

This is not a named distribution, and the normalising constant \(p(y) = \int \cdots\, d\eta\) has no closed form.

What still holds

Bayes’ rule does. The posterior exists and is perfectly well defined — we just cannot write it down.

What we do instead

Grid: evaluate on a mesh, normalise numerically
Laplace: approximate by a normal at the mode
MCMC: draw samples from it

This is the normal case

Logit regression, hierarchical models, and the horseshoe and spike-and-slab priors below are all non-conjugate.

Non-Conjugate Updating, Live

viewof nc_m = Inputs.range([-3, 3], {step: 0.1, value: -1, label: "prior mean m (log-odds)"})
viewof nc_s = Inputs.range([0.2, 3], {step: 0.05, value: 3, label: "prior sd s"})
viewof nc_n = Inputs.range([1, 200], {step: 1, value: 6, label: "observations n"})
viewof nc_p = Inputs.range([0, 1], {step: 0.01, value: 1, label: "observed share k/n"})

The pink posterior is computed on a grid of 800 points — no formula anywhere.

The green curve is the Laplace approximation: a normal centred at the posterior mode, with the curvature there. Small \(n\) with extreme shares shows its limits; large \(n\) makes it nearly exact (Bernstein–von Mises again).

The default, 6 successes in 6, has no MLE: the likelihood keeps rising, and only the prior stops the posterior — with a long right tail the normal cannot follow.

{
  const C = ldk.C;
  const n = nc_n, k = Math.round(nc_p * n), m = nc_m, s = nc_s;
  const sig = e => 1 / (1 + Math.exp(-e));
  const lp = e => k * e - n * Math.log1p(Math.exp(e)) - (e - m) * (e - m) / (2 * s * s);

  const g = ldk.seq(-7, 7, 800), dx = g[1] - g[0];
  const lv = g.map(lp), lmax = Math.max(...lv);
  let Z = 0; const un = lv.map(v => { const w = Math.exp(v - lmax); Z += w * dx; return w; });
  const post = un.map(w => w / Z);
  let mTheta = 0; for (let i = 0; i < g.length; i++) mTheta += sig(g[i]) * post[i] * dx;

  // Laplace: the mode (bisection on the score, which is decreasing), curvature there
  let lo = -30, hi = 30;
  for (let it = 0; it < 100; it++) {
    const c = (lo + hi) / 2;
    if (k - n * sig(c) - (c - m) / (s * s) > 0) lo = c; else hi = c;
  }
  const e = (lo + hi) / 2;
  const q0 = sig(e), lsd = Math.sqrt(1 / (n * q0 * (1 - q0) + 1 / (s * s)));

  const lik = g.map(v => Math.exp(k * v - n * Math.log1p(Math.exp(v))));
  const likmax = Math.max(...lik);
  const prior = g.map(v => ldk.dnorm(v, m, s));
  const ymax = 1.15 * Math.max(...post, ...prior, ldk.dnorm(0, 0, lsd));

  const { svg, p } = ldk.single(900, 500, {
    xdom: [-7, 7], ydom: [0, ymax], xlab: "log-odds η = log(θ / (1 − θ))",
    ylab: "density", margin: { t: 74 }, nx: 7
  });
  p.fill(g.map((v, i) => [v, lik[i] / likmax * ymax * 0.85]), 0, { fill: ldk.alpha(C.grey, 0.16) });
  p.path(g.map((v, i) => [v, prior[i]]), { stroke: C.blue, "stroke-width": 2.6, "stroke-dasharray": "7,5" });
  p.fill(g.map((v, i) => [v, post[i]]), 0, { fill: ldk.alpha(C.pink, 0.2) });
  p.path(g.map((v, i) => [v, post[i]]), { stroke: C.pink, "stroke-width": 3.6 });
  p.path(g.map(v => [v, ldk.dnorm(v, e, lsd)]),
         { stroke: C.green, "stroke-width": 2.6, "stroke-dasharray": "3,4" });
  p.legend([{ label: `prior  N(${ldk.fmt(m, 1)}, ${ldk.fmt(s, 2)}²) on η`, color: C.blue, lw: 2.6, dash: "7,5" },
            { label: "scaled likelihood", swatch: ldk.alpha(C.grey, 0.16) },
            { label: "posterior (grid)", color: C.pink, lw: 3.6 },
            { label: "Laplace approximation", color: C.green, lw: 2.6, dash: "3,4" }],
           { size: 14, bg: true, corner: "topleft" });
  p.header(`data ${k} of ${n}    →    posterior mean of θ = ${ldk.fmt(mTheta, 3)}`,
           { line: 1, size: 18, fill: C.pink });
  p.header(`Laplace: mode η = ${ldk.fmt(e, 2)} (θ = ${ldk.fmt(q0, 3)}), sd ${ldk.fmt(lsd, 2)}`,
           { size: 15, fill: C.green });
  return svg;
}

Shrinkage

The Posterior Mean Shrinks

Rewrite the normal–normal posterior mean as a weighted average:

\[\operatorname{E}[\theta \mid x] \;=\; w\,\mu_0 + (1-w)\,\bar x, \qquad w = \frac{1/\tau_0^2}{1/\tau_0^2 + n/\sigma^2} = \frac{\text{prior precision}}{\text{total precision}}\]

The estimate is pulled toward the prior mean by the fraction \(w\) — it is shrunk
In the strong-prior example \(w = 4/14 = 0.29\): the raw estimate 1.6 moves 29% of the way to zero, to 1.14
More data (\(n \uparrow\)) or a vaguer prior (\(\tau_0 \uparrow\)) means less shrinkage

The noisiest estimates move the most

With the same prior, a group estimated from 5 observations is pulled far toward \(\mu_0\); one estimated from 500 barely moves.

One Formula, Four Names

In estimation

Ridge regression — the Gaussian prior behind deck 04’s \(c\bar X\), shrinking coefficients toward zero
Empirical Bayes — shrinking group estimates toward the grand mean

And elsewhere

The Kalman filter — each update is exactly this, with the prior carried forward in time
Credibility theory — how insurers price a small policyholder against the portfolio

Why economists should care

Regional unemployment from 40 respondents, school effects from one cohort, firm productivity from three years — all are noisy estimates of a quantity that resembles its neighbours. Shrinkage is the answer, and it is not a fudge: it lowers mean squared error.

Shrinkage, Live

viewof nn_mu0 = Inputs.range([-3, 3], {step: 0.05, value: 0, label: "Prior mean μ₀"})
viewof nn_tau = Inputs.range([0.1, 3], {step: 0.02, value: 0.5, label: "Prior sd τ₀"})
viewof nn_xbar = Inputs.range([-3, 3], {step: 0.05, value: 1.6, label: "Sample mean x̄"})
viewof nn_n = Inputs.range([1, 200], {step: 1, value: 10, label: "n"})

The posterior always sits between the prior mean and the sample mean, and it is always narrower than both.

Shrink \(\tau_0\) toward 0.1 for a confident prior; push \(n\) up to let the data win.

{
  const C = ldk.C;
  const sigma = 1, n = nn_n, tau = nn_tau;
  const precPrior = 1 / (tau * tau), precData = n / (sigma * sigma);
  const w = precPrior / (precPrior + precData);
  const postMean = w * nn_mu0 + (1 - w) * nn_xbar;
  const postSd = Math.sqrt(1 / (precPrior + precData));
  const seData = sigma / Math.sqrt(n);

  const dom = [-4, 4];
  const g = ldk.seq(dom[0], dom[1], 600);
  const ymax = Math.max(ldk.dnorm(0, 0, postSd), ldk.dnorm(0, 0, tau),
                        ldk.dnorm(0, 0, seData));

  const { svg, p } = ldk.single(880, 520, {
    xdom: dom, ydom: [0, ymax * 1.12],
    xlab: "θ", ylab: "density", margin: { t: 74 }, nx: 5
  });

  p.path(g.map(v => [v, ldk.dnorm(v, nn_mu0, tau)]),
         { stroke: C.blue, "stroke-width": 2.6, "stroke-dasharray": "7,5" });
  p.path(g.map(v => [v, ldk.dnorm(v, nn_xbar, seData)]),
         { stroke: C.orange, "stroke-width": 2.6, "stroke-dasharray": "7,5" });
  p.fill(g.map(v => [v, ldk.dnorm(v, postMean, postSd)]), 0,
         { fill: ldk.alpha(C.pink, 0.2) });
  p.path(g.map(v => [v, ldk.dnorm(v, postMean, postSd)]),
         { stroke: C.pink, "stroke-width": 3.6 });

  p.seg(nn_mu0, 0, nn_mu0, ymax * 0.06, { stroke: C.blue, "stroke-width": 2.5 });
  p.seg(nn_xbar, 0, nn_xbar, ymax * 0.06, { stroke: C.orange, "stroke-width": 2.5 });
  p.seg(postMean, 0, postMean, ymax * 0.06, { stroke: C.pink, "stroke-width": 3 });

  p.legend([
    { label: `prior   N(${ldk.fmt(nn_mu0, 2)}, ${ldk.fmt(tau, 2)}²)`,
      color: C.blue, lw: 2.6, dash: "7,5" },
    { label: `likelihood   x̄ = ${ldk.fmt(nn_xbar, 2)}, se ${ldk.fmt(seData, 3)}`,
      color: C.orange, lw: 2.6, dash: "7,5" },
    { label: `posterior   mean ${ldk.fmt(postMean, 3)}, sd ${ldk.fmt(postSd, 3)}`,
      color: C.pink, lw: 3.6 }
  ], { size: 14, bg: true });

  p.header(`weight on the prior  w = ${ldk.fmt(w, 3)}` +
           `      weight on the data  ${ldk.fmt(1 - w, 3)}`,
           { line: 1, size: 18, fill: C.pink });
  p.header(`posterior precision ${ldk.fmt(precPrior + precData, 2)} = ` +
           `${ldk.fmt(precPrior, 2)} (prior) + ${ldk.fmt(precData, 2)} (data)`,
           { size: 15, fill: C.grey });
  return svg;
}

Stronger Shrinkage: The Horseshoe Prior

A normal prior shrinks every estimate by the same fraction — noise and large effects alike. With many coefficients, most near zero, we want to shrink the small ones hard and leave the large ones alone.

\[\theta_j \mid \lambda_j, \tau \;\sim\; N(0,\ \lambda_j^2\,\tau^2)\] \[\lambda_j \sim \text{Cauchy}^{+}(0, 1)\]

\(\tau\) — global: how much shrinkage overall
\(\lambda_j\) — local: one per coefficient, heavy-tailed, so a single \(\theta_j\) can escape

Why “horseshoe”

The shrinkage factor \(\kappa_j = 1/(1+\lambda_j^2\tau^2)\) piles up at 0 (keep) and 1 (kill): a horseshoe-shaped density.

{
  const C = ldk.C;
  // horseshoe density with tau = 1: substitute lambda = tan(pi u / 2), u ~ U(0,1)
  const U = ldk.seq(0.00025, 0.99975, 2000);
  const hs = t => { let s = 0; for (const u of U) s += ldk.dnorm(t, 0, Math.tan(Math.PI * u / 2)); return s / U.length; };
  const lap = t => Math.exp(-Math.abs(t) * Math.SQRT2) / Math.SQRT2;      // variance 1
  const svg = ldk.chart(760, 480);

  const g1 = ldk.seq(-3, 3, 241);
  const p1 = ldk.panel(svg, { left: 0, top: 0, w: 760, h: 240, xdom: [-3, 3], ydom: [0, 1.2],
    xlab: "θ", ylab: "prior density", margin: { l: 60, r: 14, t: 36, b: 46 }, nx: 6, ny: 3 });
  p1.path(g1.map(t => [t, ldk.dnorm(t, 0, 1)]), { stroke: C.blue, "stroke-width": 3 });
  p1.path(g1.map(t => [t, lap(t)]), { stroke: C.orange, "stroke-width": 3, "stroke-dasharray": "7,5" });
  p1.path(g1.map(t => [t, hs(t)]), { stroke: C.pink, "stroke-width": 3.4 });
  p1.legend([{ label: "normal (ridge)", color: C.blue, lw: 3 },
             { label: "Laplace (lasso)", color: C.orange, lw: 3, dash: "7,5" },
             { label: "horseshoe", color: C.pink, lw: 3.4 }], { size: 14, bg: true });
  p1.header("near zero: the horseshoe has an infinitely tall spike", { size: 15 });

  const g2 = ldk.seq(2, 10, 161);
  const p2 = ldk.panel(svg, { left: 0, top: 240, w: 380, h: 240, xdom: [2, 10], ydom: [-7, -0.5],
    xlab: "θ", ylab: "log₁₀ density", margin: { l: 60, r: 14, t: 36, b: 46 }, nx: 4, ny: 3 });
  p2.path(g2.map(t => [t, Math.log10(ldk.dnorm(t, 0, 1))]), { stroke: C.blue, "stroke-width": 3 });
  p2.path(g2.map(t => [t, Math.log10(lap(t))]), { stroke: C.orange, "stroke-width": 3, "stroke-dasharray": "7,5" });
  p2.path(g2.map(t => [t, Math.log10(hs(t))]), { stroke: C.pink, "stroke-width": 3.4 });
  p2.header("far out: heavy tails", { size: 15 });

  const g3 = ldk.seq(0.01, 0.99, 197);
  const p3 = ldk.panel(svg, { left: 380, top: 240, w: 380, h: 240, xdom: [0, 1], ydom: [0, 3.5],
    xlab: "shrinkage factor κ", margin: { l: 40, r: 14, t: 36, b: 46 }, nx: 4, yticks: [] });
  p3.fill(g3.map(k => [k, ldk.dbeta(k, 0.5, 0.5)]), 0, { fill: ldk.alpha(C.pink, 0.2) });
  p3.path(g3.map(k => [k, ldk.dbeta(k, 0.5, 0.5)]), { stroke: C.pink, "stroke-width": 3.4 });
  p3.vline(0.5, { stroke: C.blue, "stroke-width": 3 });
  p3.header("the horseshoe: Beta(½, ½)", { size: 15, fill: C.pink });
  p3.add(ldk.el("text", { x: p3.x(0.03), y: p3.y(3.1), "font-size": 13, fill: C.grey }, "keep"), true);
  p3.add(ldk.el("text", { x: p3.x(0.97), y: p3.y(3.1), "font-size": 13, fill: C.grey,
                          "text-anchor": "end" }, "shrink to 0"), true);
  p3.add(ldk.el("text", { x: p3.x(0.52), y: p3.y(1.6), "font-size": 13, fill: C.blue }, "normal prior: one fixed κ"), true);
  return svg;
}

Horseshoe versus Normal, Live

viewof hs_y = Inputs.range([-8, 8], {step: 0.1, value: 1.2, label: "observed estimate y"})
viewof hs_tau = Inputs.range([0.05, 2], {step: 0.05, value: 0.3, label: "horseshoe global τ"})
viewof hs_sd = Inputs.range([0.2, 3], {step: 0.05, value: 1, label: "normal prior sd"})

One noisy estimate \(y \sim N(\theta, 1)\) — think of one coefficient among many.

Left: the posterior mean as a function of \(y\). The normal prior shrinks by the same fraction everywhere: a straight line. The horseshoe crushes small \(y\) to about zero but follows the diagonal once \(y\) is clearly large.

Right: the two posteriors for the current \(y\). Move \(y\) from 1 to 6.

{
  const C = ldk.C;
  const tau = hs_tau, s = hs_sd, y = hs_y;
  const U = ldk.seq(0.0005, 0.9995, 1000);
  const lam2 = U.map(u => { const l = Math.tan(Math.PI * u / 2); return l * l * tau * tau; });

  // horseshoe posterior mean: y * E[1 - kappa | y]
  const hsMean = yy => {
    let num = 0, den = 0;
    for (const v of lam2) { const w = ldk.dnorm(yy, 0, Math.sqrt(1 + v)); num += w * v / (1 + v); den += w; }
    return yy * num / den;
  };
  const shr = s * s / (1 + s * s);
  const ys = ldk.seq(-8, 8, 161);

  const svg = ldk.chart(920, 480);
  const pl = ldk.panel(svg, { left: 0, w: 460, xdom: [-8, 8], ydom: [-8, 8],
    xlab: "observed y", ylab: "posterior mean E[θ | y]", margin: { l: 62, r: 14, t: 74, b: 50 }, nx: 4, ny: 4 });
  pl.path([[-8, -8], [8, 8]], { stroke: ldk.alpha(C.grey, 0.6), "stroke-width": 1.6, "stroke-dasharray": "5,5" });
  pl.path(ys.map(v => [v, v * shr]), { stroke: C.blue, "stroke-width": 3 });
  pl.path(ys.map(v => [v, hsMean(v)]), { stroke: C.pink, "stroke-width": 3.4 });
  const mh = hsMean(y), mn = y * shr;
  pl.dot(y, mn, 7, { fill: C.blue, stroke: "white", "stroke-width": 2 });
  pl.dot(y, mh, 7, { fill: C.pink, stroke: "white", "stroke-width": 2 });
  pl.legend([{ label: "no shrinkage", color: C.grey, lw: 1.6, dash: "5,5" },
             { label: "normal prior", color: C.blue, lw: 3 },
             { label: "horseshoe", color: C.pink, lw: 3.4 }], { size: 13, bg: true, corner: "topleft" });
  pl.header(`y = ${ldk.fmt(y, 1)}:  normal → ${ldk.fmt(mn, 2)},  horseshoe → ${ldk.fmt(mh, 2)}`,
            { line: 1, size: 17 });
  const kept = v => Math.abs(y) < 0.05 ? "—" : `${ldk.fmt(100 * v / y, 0)}%`;
  pl.header(`share of y kept:  normal ${kept(mn)},  horseshoe ${kept(mh)}`, { size: 14, fill: C.grey });

  // right: posterior densities at the current y
  const g = ldk.seq(-9, 9, 721);
  const hsPrior = g.map(t => { let a = 0; for (const v of lam2) a += ldk.dnorm(t, 0, Math.sqrt(v)); return a / lam2.length; });
  const dx = g[1] - g[0];
  let Z = 0; const ph = g.map((t, i) => { const w = ldk.dnorm(y, t, 1) * hsPrior[i]; Z += w * dx; return w; });
  const postH = ph.map(w => w / Z);
  const pnSd = Math.sqrt(shr), postN = g.map(t => ldk.dnorm(t, mn, pnSd));
  const ymax = 1.15 * Math.max(...postN, 0.6);
  const pr = ldk.panel(svg, { left: 460, w: 460, xdom: [-9, 9], ydom: [0, ymax],
    xlab: "θ", ylab: "posterior density", margin: { l: 58, r: 14, t: 74, b: 50 }, nx: 6, yticks: [] });
  pr.vline(y, { stroke: ldk.alpha(C.grey, 0.6), "stroke-width": 1.6, "stroke-dasharray": "5,5" });
  pr.fill(g.map((t, i) => [t, postN[i]]), 0, { fill: ldk.alpha(C.blue, 0.15) });
  pr.path(g.map((t, i) => [t, postN[i]]), { stroke: C.blue, "stroke-width": 3 });
  pr.fill(g.map((t, i) => [t, postH[i]]), 0, { fill: ldk.alpha(C.pink, 0.18) });
  pr.path(g.map((t, i) => [t, postH[i]]), { stroke: C.pink, "stroke-width": 3.4 });
  pr.header("posteriors at this y", { line: 1, size: 17 });
  pr.header("dashed: the observed y", { size: 14, fill: C.grey });
  return svg;
}

When to Reach for the Horseshoe

Sparse problems

Many candidate regressors, few that matter: growth determinants, many controls
Many group effects, most negligible: regional or firm effects
Forecasting with large VARs

In practice

Non-conjugate — fitted by MCMC
R: brms with prior = horseshoe(), or rstanarm with prior = hs()
Its posterior mean is never exactly zero — to select variables you need the next idea

Ridge, lasso, horseshoe

All three are priors: ridge is the normal, lasso the Laplace (at its mode), and the horseshoe gives selective shrinkage — none for large effects, almost total for small ones.

Two Kinds of Interval

Credible versus Confidence

Credible interval

A set \(C\) with \[P(\theta \in C \mid y) = 1-\alpha\]

*“Given the data and the prior, there is a 95% probability that* \(\theta\) lies in here.”

Confidence interval

A rule with \[P_\theta\big(\theta \in C(Y)\big) = 1-\alpha \ \ \forall\theta\]

*“95% of intervals built this way contain* \(\theta\).”

The credible interval is what people think a confidence interval is

Deck 05 listed \(P(H_0 \mid \text{data})\) as something a p-value is not, and deck 04 listed the probability statement as something a confidence interval is not. The Bayesian posterior is the object that actually licenses it — at the cost of requiring a prior.

Two Intervals, One Dataset

viewof cc_n = Inputs.range([5, 200], {step: 1, value: 20, label: "n"})
viewof cc_k = Inputs.range([0, 40], {step: 1, value: 1, label: "successes k"})
viewof cc_prior = Inputs.select(
  new Map([["Uniform — Beta(1,1)", "flat"],
           ["Jeffreys — Beta(½,½)", "jeff"],
           ["Sceptical — Beta(2,18)", "scept"]]),
  {label: "Prior", value: "flat"})

With a flat prior and a large \(n\) the two intervals nearly coincide — which is why the misinterpretation survives.

Push \(k\) toward 0 or \(n\), or shrink \(n\), and they separate. Only one of them stays inside \([0,1]\) by construction.

{
  const C = ldk.C;
  const n = cc_n, k = Math.min(cc_k, n);
  const pri = { flat: [1, 1], jeff: [0.5, 0.5], scept: [2, 18] }[cc_prior];
  const a = pri[0] + k, b = pri[1] + n - k;

  const credLo = ldk.qbeta(0.025, a, b), credHi = ldk.qbeta(0.975, a, b);
  const phat = k / n, z = 1.959964;
  const s = Math.sqrt(phat * (1 - phat) / n);
  const waldLo = phat - z * s, waldHi = phat + z * s;
  const d = 1 + z * z / n;
  const wc = (phat + z * z / (2 * n)) / d;
  const wh = z * Math.sqrt(phat * (1 - phat) / n + z * z / (4 * n * n)) / d;

  const hi = Math.min(1, Math.max(credHi, waldHi, wc + wh, phat) + 0.10);
  const lo = Math.max(-0.12, Math.min(credLo, waldLo, wc - wh, phat) - 0.06);

  const g = ldk.seq(Math.max(lo, 1e-4), Math.min(hi, 0.9999), 500);
  const dens = g.map(t => ldk.dbeta(t, a, b));
  const ymax = Math.max(...dens);

  const { svg, p } = ldk.single(880, 520, {
    xdom: [lo, hi], ydom: [-ymax * 0.30, ymax * 1.10],
    xlab: "θ", ylab: "posterior density", margin: { t: 74 }, nx: 6,
    yticks: []
  });

  const inside = g.map((t, i) => [t, t >= credLo && t <= credHi ? dens[i] : 0]);
  p.fill(inside, 0, { fill: ldk.alpha(C.pink, 0.28) });
  p.path(g.map((t, i) => [t, dens[i]]), { stroke: C.pink, "stroke-width": 3.4 });

  const bar = (y, l, h, col, lab) => {
    p.seg(l, y, h, y, { stroke: col, "stroke-width": 6 });
    p.seg(l, y - ymax * 0.022, l, y + ymax * 0.022, { stroke: col, "stroke-width": 2.5 });
    p.seg(h, y - ymax * 0.022, h, y + ymax * 0.022, { stroke: col, "stroke-width": 2.5 });
    p.add(ldk.el("text", { x: p.x(lo) + 6, y: p.y(y) - 6, "font-size": 14,
                           fill: col }, lab), true);
  };
  bar(-ymax * 0.09, credLo, credHi, C.pink, "credible (Bayes)");
  bar(-ymax * 0.18, wc - wh, wc + wh, C.blue, "Wilson");
  bar(-ymax * 0.27, waldLo, waldHi, C.orange, "Wald");
  p.vline(0, { stroke: ldk.alpha(C.grey, 0.55), "stroke-width": 1.6,
               "stroke-dasharray": "4,5" });

  p.header(`credible [${ldk.fmt(credLo, 3)}, ${ldk.fmt(credHi, 3)}]` +
           `      Wald [${ldk.fmt(waldLo, 3)}, ${ldk.fmt(waldHi, 3)}]`,
           { line: 1, size: 17, fill: C.pink });
  p.header(waldLo < 0 ? "the Wald interval has left the parameter space"
                      : `posterior mean ${ldk.fmt(a / (a + b), 3)},  k/n = ${ldk.fmt(phat, 3)}`,
           { size: 15, fill: waldLo < 0 ? C.pink : C.grey });
  return svg;
}

The Two Sentences

With a flat prior and plenty of data the intervals nearly coincide numerically — but they are answers to different questions
A credible interval is a statement about \(\theta\) given these data. A confidence interval is a statement about the procedure across hypothetical repetitions
The posterior respects the parameter space: a probability cannot be negative, and the credible interval never is

The honest trade

The Bayesian statement is the one you wanted. It costs you a prior, and someone can always ask where the prior came from.

Evidence

Bayes Factors

Definition — Bayes Factor

\[\text{BF}_{10} \;=\; \frac{p(y \mid H_1)}{p(y \mid H_0)} \;=\; \frac{\int p(y\mid\theta)\,p_1(\theta)\,d\theta}{p(y\mid\theta_0)}\] It converts prior odds into posterior odds: \[\frac{P(H_1\mid y)}{P(H_0\mid y)} \;=\; \text{BF}_{10}\times\frac{P(H_1)}{P(H_0)}\]

What it does that a p-value cannot

It can favour \(H_0\). A p-value of 0.6 is not evidence for the null; a Bayes factor of \(\tfrac{1}{12}\) is.

What it costs

The integral needs a prior under \(H_1\) — and unlike the posterior, the Bayes factor never stops depending on it.

Model Selection: Spike and Slab

Many Models at Once

A Bayes factor compares two models. With \(p\) candidate regressors there are \(2^p\) models: 1,024 for \(p = 10\), a billion for \(p = 30\).

Give every coefficient an indicator \(\gamma_j \in \{0, 1\}\): in the model or out
Put a prior on the indicators, and let the posterior decide which regressors belong
Read off each regressor’s posterior inclusion probability

\[\text{PIP}_j \;=\; P(\gamma_j = 1 \mid y)\]

What you report

The PIPs instead of significance stars, and the median probability model: every regressor with \(\text{PIP}_j > 0.5\).

Where economists use it

Growth regressions with dozens of candidate determinants, and forecasting with many predictors: Bayesian model averaging.

The Spike-and-Slab Prior

\[\beta_j \mid \gamma_j \;\sim\; (1-\gamma_j)\,\underbrace{\delta_0}_{\text{spike}} \;+\; \gamma_j\,\underbrace{N(0, c_1^2)}_{\text{slab}}\] \[\gamma_j \;\sim\; \text{Bernoulli}(\pi)\]

Spike \(\delta_0\): a point mass — the coefficient is exactly zero
Slab (\(c_1\) wide): “some plausible non-zero value”
\(\pi\): the prior share of regressors you expect to matter

SSVS — stochastic search variable selection

A Gibbs sampler alternates \(\beta \mid \gamma\) (a conjugate normal update for the included coefficients; the rest are exactly 0) and \(\gamma_j \mid y, \gamma_{-j}\) (a coin flip per regressor, with \(\beta_j\) integrated out). \(\text{PIP}_j\) = share of draws with \(\gamma_j = 1\).

{
  const C = ldk.C;
  const c1 = 1.5, pi = 0.5;
  const g = ldk.seq(-4.5, 4.5, 901);
  const { svg, p } = ldk.single(700, 500, {
    xdom: [-4.5, 4.5], ydom: [0, 0.6], xlab: "β",
    ylab: "prior density  /  point mass", margin: { t: 64 }, nx: 6, ny: 3
  });
  const slab = g.map(b => pi * ldk.dnorm(b, 0, c1));
  p.fill(g.map((b, i) => [b, slab[i]]), 0, { fill: ldk.alpha(C.orange, 0.25) });
  p.path(g.map((b, i) => [b, slab[i]]), { stroke: C.orange, "stroke-width": 3 });
  p.spike(0, 1 - pi, { color: C.blue, lw: 4 });
  p.add(ldk.el("text", { x: p.x(0.18), y: p.y(1 - pi) + 6, "font-size": 15, fill: C.blue },
               `mass 1 − π = ${1 - pi} at exactly 0`), true);
  p.legend([{ label: "spike  (1 − π) · δ₀", color: C.blue, lw: 4 },
            { label: `slab  π · N(0, ${c1}²)`, color: C.orange, lw: 3 }],
           { size: 15, bg: true });
  p.header(`π = ${pi}: half the regressors are expected to matter`, { line: 1, size: 17 });
  p.header("arrow height = probability, curve = density", { size: 14, fill: C.grey });
  return svg;
}

Inclusion Probability, Live

One coefficient, estimated as \(\hat\beta\) with standard error 1. Is it spike or slab?

viewof ss_t = Inputs.range([0, 6], {step: 0.05, value: 2, label: "estimate (= t-statistic)"})
viewof ss_c1 = Inputs.range([0.5, 10], {step: 0.1, value: 3, label: "slab sd c₁"})
viewof ss_pi = Inputs.range([0.05, 0.95], {step: 0.01, value: 0.5, label: "prior inclusion π"})

The PIP is prior odds times a Bayes factor — slab against spike — turned back into a probability:

\(\dfrac{\text{PIP}}{1-\text{PIP}} = \dfrac{\pi}{1-\pi} \cdot \dfrac{N(\hat\beta;\,0,\,1+c_1^2)}{N(\hat\beta;\,0,\,1)}\)

So the Bayes-factor lesson carries over: the answer depends on the slab width \(c_1\). Make the slab very wide and a \(t\) of 2 no longer earns inclusion.

{
  const C = ldk.C;
  const t = ss_t, c1 = ss_c1, pi = ss_pi;
  const pip = x => {
    const a = pi * ldk.dnorm(x, 0, Math.sqrt(1 + c1 * c1));
    const b = (1 - pi) * ldk.dnorm(x, 0, 1);
    return a / (a + b);
  };
  const P = pip(t);
  const bf = ldk.dnorm(t, 0, Math.sqrt(1 + c1 * c1)) / ldk.dnorm(t, 0, 1);

  const svg = ldk.chart(920, 480);
  // left: prior and posterior for beta, each a point mass at 0 plus a normal slab
  const g = ldk.seq(-3, 7, 801);
  const v = c1 * c1 / (1 + c1 * c1);
  const prior = g.map(b => pi * ldk.dnorm(b, 0, c1));
  const pSlab = g.map(b => P * ldk.dnorm(b, t * v, Math.sqrt(v)));
  const pl = ldk.panel(svg, { left: 0, w: 480, xdom: [-3, 7], ydom: [0, 1.2],
    xlab: "β", ylab: "density", margin: { l: 58, r: 14, t: 74, b: 50 }, nx: 5, ny: 3 });
  pl.fill(g.map((b, i) => [b, prior[i]]), 0, { fill: ldk.alpha(C.grey, 0.14) });
  pl.fill(g.map((b, i) => [b, pSlab[i]]), 0, { fill: ldk.alpha(C.orange, 0.25) });
  pl.path(g.map((b, i) => [b, pSlab[i]]), { stroke: C.orange, "stroke-width": 3 });
  pl.path(g.map((b, i) => [b, prior[i]]), { stroke: C.grey, "stroke-width": 2.6, "stroke-dasharray": "6,5" });
  pl.spike(-0.1, 1 - pi, { color: C.grey, lw: 2.6, dash: "6,5", head: 6 });
  pl.spike(0.1, 1 - P, { color: C.blue, lw: 4 });
  pl.vline(t, { stroke: C.pink, "stroke-width": 2, "stroke-dasharray": "4,4" });
  pl.legend([{ label: `prior: slab π · N(0, ${c1}²) and spike`, color: C.grey, lw: 2.6, dash: "6,5" },
             { label: `posterior: slab part (${ldk.fmt(P, 2)})`, color: C.orange, lw: 3 },
             { label: `posterior: spike δ₀ (mass ${ldk.fmt(1 - P, 2)})`, color: C.blue, lw: 4 }],
            { size: 13, bg: true });
  pl.header(`PIP = ${ldk.fmt(P, 3)}`, { line: 1, size: 19, fill: P > 0.5 ? C.green : C.pink });
  pl.header(`Bayes factor slab : spike = ${bf >= 100 ? bf.toExponential(1) : ldk.fmt(bf, 2)}`,
            { size: 14, fill: C.grey });

  // right: PIP as a function of the t-statistic
  const ts = ldk.seq(0, 6, 241);
  const pr = ldk.panel(svg, { left: 480, w: 440, xdom: [0, 6], ydom: [0, 1],
    xlab: "|β̂| / se   (the t-statistic)", ylab: "PIP", margin: { l: 58, r: 14, t: 74, b: 50 }, nx: 6, ny: 4 });
  pr.hline(0.5, { stroke: ldk.alpha(C.grey, 0.7), "stroke-width": 1.6, "stroke-dasharray": "5,5" });
  pr.vline(1.96, { stroke: ldk.alpha(C.orange, 0.8), "stroke-width": 1.6, "stroke-dasharray": "5,5" });
  pr.path(ts.map(x => [x, pip(x)]), { stroke: C.ink, "stroke-width": 3 });
  pr.dot(t, P, 8, { fill: C.pink, stroke: "white", "stroke-width": 2 });
  pr.add(ldk.el("text", { x: pr.x(2.02), y: pr.y(0.06), "font-size": 13, fill: C.orange }, "p = 0.05"), true);
  pr.header(P > 0.5 ? "in the median probability model" : "left out of the median probability model",
            { line: 1, size: 17, fill: P > 0.5 ? C.green : C.pink });
  pr.header("dashed: PIP = ½ and |t| = 1.96", { size: 14, fill: C.grey });
  return svg;
}

Choosing Regressors, Live

Eight candidate regressors; only \(x_1, x_2, x_3\) matter (\(\beta = 0.6,\ 0.35,\ 0.2\)), and \(x_4\) is correlated with \(x_1\) (\(\rho = 0.7\)).

viewof sv_n = Inputs.range([20, 500], {step: 10, value: 80, label: "sample size n"})
viewof sv_pi = Inputs.range([0.05, 0.95], {step: 0.05, value: 0.5, label: "prior inclusion π"})
viewof sv_seed = Inputs.range([1, 30], {step: 1, value: 5, label: "data set #"})

With a point-mass spike and a \(g\)-prior slab (\(g = n\)), all \(2^8 = 256\) models can be scored exactly in the browser. SSVS samples the same posterior when \(2^p\) is too large to enumerate.

Bars: \(\text{PIP}_j\). A ★ marks \(|t| > 1.96\) in the full OLS regression.

Try small \(n\): the weak \(x_3\) drifts in and out, and \(x_4\) borrows credit from \(x_1\).

{
  const C = ldk.C;
  const p = 8, n = sv_n, pi = sv_pi, g = n;
  const beta = [0.6, 0.35, 0.2, 0, 0, 0, 0, 0];
  const u = ldk.rng(1000 + sv_seed);

  // simulate
  const X = [], y = [];
  for (let i = 0; i < n; i++) {
    const r = [];
    for (let j = 0; j < p; j++) r.push(u.normal());
    r[3] = 0.7 * r[0] + Math.sqrt(1 - 0.49) * r[3];
    X.push(r);
    y.push(r.reduce((s, v, j) => s + v * beta[j], 0) + u.normal());
  }
  // centre
  const mx = ldk.seq(0, p - 1, p).map(j => ldk.mean(X.map(r => r[j])));
  const my = ldk.mean(y);
  const Xc = X.map(r => r.map((v, j) => v - mx[j])), yc = y.map(v => v - my);
  const XtX = [], Xty = [];
  for (let a = 0; a < p; a++) {
    XtX.push([]); let s = 0;
    for (let i = 0; i < n; i++) s += Xc[i][a] * yc[i];
    Xty.push(s);
    for (let b = 0; b < p; b++) { let q = 0; for (let i = 0; i < n; i++) q += Xc[i][a] * Xc[i][b]; XtX[a].push(q); }
  }
  const yty = yc.reduce((s, v) => s + v * v, 0);
  const solve = (A, b) => {                     // Gaussian elimination, small systems
    const m = A.length, M = A.map((r, i) => [...r, b[i]]);
    for (let c = 0; c < m; c++) {
      let piv = c; for (let r = c + 1; r < m; r++) if (Math.abs(M[r][c]) > Math.abs(M[piv][c])) piv = r;
      [M[c], M[piv]] = [M[piv], M[c]];
      for (let r = c + 1; r < m; r++) { const f = M[r][c] / M[c][c]; for (let k = c; k <= m; k++) M[r][k] -= f * M[c][k]; }
    }
    const x = new Array(m).fill(0);
    for (let r = m - 1; r >= 0; r--) { let s = M[r][m]; for (let k = r + 1; k < m; k++) s -= M[r][k] * x[k]; x[r] = s / M[r][r]; }
    return x;
  };

  // score every model
  const models = [];
  for (let mask = 0; mask < (1 << p); mask++) {
    const idx = []; for (let j = 0; j < p; j++) if (mask & (1 << j)) idx.push(j);
    let R2 = 0;
    if (idx.length) {
      const b = solve(idx.map(a => idx.map(c => XtX[a][c])), idx.map(a => Xty[a]));
      R2 = b.reduce((s, v, k) => s + v * Xty[idx[k]], 0) / yty;
    }
    const k = idx.length;
    const lml = 0.5 * (n - 1 - k) * Math.log(1 + g) - 0.5 * (n - 1) * Math.log(1 + g * (1 - R2));
    models.push({ mask, idx, lp: lml + k * Math.log(pi) + (p - k) * Math.log(1 - pi) });
  }
  const lmax = Math.max(...models.map(m => m.lp));
  let Z = 0; for (const m of models) { m.w = Math.exp(m.lp - lmax); Z += m.w; }
  for (const m of models) m.w /= Z;
  const pip = ldk.seq(0, p - 1, p).map(j => models.reduce((s, m) => s + (m.mask & (1 << j) ? m.w : 0), 0));

  // full OLS t-statistics
  const bFull = solve(XtX, Xty);
  const s2 = (yty - bFull.reduce((s, v, j) => s + v * Xty[j], 0)) / (n - 1 - p);
  const tstat = ldk.seq(0, p - 1, p).map(j => {
    const e = new Array(p).fill(0); e[j] = 1;
    return bFull[j] / Math.sqrt(s2 * solve(XtX, e)[j]);
  });

  const svg = ldk.chart(940, 500);
  const pl = ldk.panel(svg, { left: 0, w: 600, xdom: [0.4, 8.6], ydom: [0, 1.12],
    xlab: "regressor", ylab: "posterior inclusion probability", margin: { l: 62, r: 14, t: 70, b: 50 },
    xticks: [1, 2, 3, 4, 5, 6, 7, 8], xfmt: v => `x${v}`, yticks: [0, 0.25, 0.5, 0.75, 1] });
  pl.hline(0.5, { stroke: ldk.alpha(C.grey, 0.7), "stroke-width": 1.6, "stroke-dasharray": "5,5" });
  for (let j = 0; j < p; j++) {
    pl.bar(j + 1 - 0.32, j + 1 + 0.32, 0, pip[j],
           { fill: beta[j] !== 0 ? C.green : (j === 3 ? C.orange : ldk.alpha(C.grey, 0.7)) });
    if (Math.abs(tstat[j]) > 1.96)
      pl.add(ldk.el("text", { x: pl.x(j + 1), y: pl.y(Math.min(pip[j], 1)) - 6, "font-size": 20,
                              "text-anchor": "middle", fill: C.pink }, "★"), true);
  }
  pl.legend([{ label: "true β ≠ 0", swatch: C.green },
             { label: "β = 0, correlated with x1", swatch: C.orange },
             { label: "β = 0", swatch: ldk.alpha(C.grey, 0.7) }], { size: 13, bg: true, corner: "topright" });
  const mpm = pip.map((v, j) => v > 0.5 ? `x${j + 1}` : null).filter(Boolean);
  pl.header(`median probability model:  { ${mpm.join(", ") || "none"} }`, { line: 1, size: 17, fill: C.pink });
  pl.header(`n = ${n}, π = ${ldk.fmt(pi, 2)}`, { size: 14, fill: C.grey });

  // right: the top models
  const top = [...models].sort((a, b) => b.w - a.w).slice(0, 8);
  const x0 = 620;
  svg.append(ldk.el("text", { x: x0, y: 46, "font-size": 17, "font-weight": "bold", fill: C.ink }, "most probable models"));
  top.forEach((m, i) => {
    const yy = 84 + i * 44;
    svg.append(ldk.el("rect", { x: x0, y: yy - 16, width: 160 * m.w / top[0].w, height: 22,
                                fill: ldk.alpha(C.pink, 0.25) }));
    svg.append(ldk.el("text", { x: x0 + 6, y: yy, "font-size": 15, fill: C.ink },
      `{ ${m.idx.map(j => "x" + (j + 1)).join(", ") || "none"} }`));
    svg.append(ldk.el("text", { x: 930, y: yy, "font-size": 15, fill: C.grey, "text-anchor": "end" },
      ldk.fmt(m.w, 3)));
  });
  return svg;
}

What Spike-and-Slab Buys You

Model uncertainty, stated. Not “the model”, but how much each regressor is worth including, given everything else that might be included
Averaging instead of picking. Predictions weighted by posterior model probability — Bayesian model averaging — beat any single selected model on average
Correlated regressors share the credit. When \(x_1\) and \(x_4\) move together, their PIPs trade off rather than both looking significant

Priors still matter

Like every Bayes factor, the PIPs depend on the slab width and on \(\pi\). Report them, and show what happens when they change.

In R

BoomSpikeSlab::lm.spike() (SSVS-style sampler), BAS::bas.lm() and BMS::bms() (model averaging over the model space).

Decisions

Loss, Risk, and What “Best” Means

Deck 04 defined “best” with a loss, averaged over the data. The Bayes estimator minimises the same losses, averaged over the posterior instead.

Definition — Bayes Estimator

Given a loss \(L(\theta, d)\) for reporting \(d\) when the truth is \(\theta\), the Bayes estimator minimises posterior expected loss: \[\hat\theta \;=\; \arg\min_{d}\ \operatorname{E}\big[L(\theta,d) \mid y\big] \;=\; \arg\min_d \int L(\theta,d)\,p(\theta\mid y)\,d\theta\]

Loss function	Bayes estimator
Squared error \((\theta-d)^2\)	posterior mean
Absolute error \(\lvert\theta-d\rvert\)	posterior median
0–1 loss	posterior mode (the MAP)
Asymmetric linear, costs \(a\) under / \(b\) over	posterior quantile at \(a/(a+b)\)

Asymmetric Loss, Live

Under-forecasting a deficit and over-forecasting it do not cost the same.

viewof al_ratio = Inputs.range([0.1, 10], {step: 0.1, value: 4,
                                           label: "Cost of under ÷ over"})
viewof al_skew = Inputs.range([0.2, 1.4], {step: 0.02, value: 0.8,
                                           label: "Posterior skew"})

With symmetric loss the optimum is the median. Make under-estimating four times as costly and the optimal report moves into the upper tail — deliberately, not as a bias.

The mean, median and mode of the posterior have not moved at all.

{
  const C = ldk.C;
  const s = al_skew, mu = 0;                       // lognormal posterior
  const dens = x => x <= 0 ? 0 : ldk.dlnorm(x, mu, s);
  const q = pr => Math.exp(mu + s * ldk.qnorm(pr));

  const a = al_ratio, b = 1;                       // under costs a, over costs b
  const target = a / (a + b);
  const dOpt = q(target);
  const median = q(0.5);
  const mean = Math.exp(mu + s * s / 2);

  const hi = q(0.985), g = ldk.seq(1e-3, hi, 600);
  const ys = g.map(dens), ymax = Math.max(...ys);

  // Posterior expected loss, by quadrature. The integration grid runs far into
  // the upper tail even though the plot stops at q(0.985): with asymmetric loss
  // that tail is exactly the region being penalised, and truncating it biases
  // the minimiser downwards.
  const gi = ldk.seq(1e-4, q(0.99995), 3000);
  const yi = gi.map(dens);
  const risk = d => {
    let r = 0;
    for (let i = 1; i < gi.length; i++) {
      const t = gi[i], w = (gi[i] - gi[i - 1]) * yi[i];
      r += w * (t > d ? a * (t - d) : b * (d - t));
    }
    return r;
  };
  const acts = ldk.seq(q(0.02), q(0.98), 160);
  const rs = acts.map(risk);
  const rmin = Math.min(...rs), rmax = Math.max(...rs);

  const svg = ldk.chart(880, 520);

  const pl = ldk.panel(svg, {
    left: 0, w: 880, h: 300, xdom: [0, hi], ydom: [0, ymax * 1.14],
    ylab: "posterior", margin: { l: 66, r: 20, t: 66, b: 8 },
    nx: 6, xticks: []
  });
  const lower = g.map((x, i) => [x, x <= dOpt ? ys[i] : 0]);
  pl.fill(g.map((x, i) => [x, ys[i]]), 0, { fill: ldk.alpha(C.grey, 0.13) });
  pl.fill(lower, 0, { fill: ldk.alpha(C.blue, 0.22) });
  pl.path(g.map((x, i) => [x, ys[i]]), { stroke: C.blue, "stroke-width": 3.2 });
  pl.vline(median, { stroke: C.grey, "stroke-width": 2.4, "stroke-dasharray": "6,5" });
  pl.vline(mean, { stroke: C.green, "stroke-width": 2.4, "stroke-dasharray": "6,5" });
  pl.vline(dOpt, { stroke: C.pink, "stroke-width": 4 });
  pl.legend([{ label: `posterior median ${ldk.fmt(median, 2)}`, color: C.grey, lw: 2.4, dash: "6,5" },
             { label: `posterior mean ${ldk.fmt(mean, 2)}`, color: C.green, lw: 2.4, dash: "6,5" },
             { label: `optimal report ${ldk.fmt(dOpt, 2)}`, color: C.pink, lw: 4 }],
            { size: 14, bg: true });
  pl.header(`under-estimating costs ${ldk.fmt(a, 1)}× as much as over-estimating` +
            `  →  report the ${ldk.fmt(100 * target, 0)}ᵗʰ percentile`,
            { size: 18, fill: C.pink });

  const pr2 = ldk.panel(svg, {
    left: 0, top: 300, w: 880, h: 220, xdom: [0, hi],
    ydom: [rmin - (rmax - rmin) * 0.12, rmax],
    xlab: "reported value  d", ylab: "expected loss",
    margin: { l: 66, r: 20, t: 12, b: 52 }, nx: 6, ny: 3
  });
  pr2.path(acts.map((d, i) => [d, rs[i]]), { stroke: C.orange, "stroke-width": 3.2 });
  pr2.vline(dOpt, { stroke: C.pink, "stroke-width": 3 });
  pr2.dot(dOpt, risk(dOpt), 7, { fill: C.pink, stroke: "white", "stroke-width": 2 });
  return svg;
}

Why Economists Should Recognise This

It is expected-utility maximisation, with the unknown parameter as the state of the world — the framework is already familiar from micro
The estimator is not a property of the data alone. It depends on what a mistake costs
Inventory, capacity, deficit forecasts, capital buffers and reserve prices are all asymmetric-loss problems, and reporting the posterior mean is simply the wrong answer for them

The unifying claim

“Best” is never a purely statistical question. Squared-error loss is a choice, and it is the one that happens to be conventional.

Computation and Practice

The Computation Ladder

Conjugate — closed form, as in the Beta–Bernoulli and normal–normal models. Fast, exact, and restricted to a handful of model families
Grid approximation — evaluate the posterior on a mesh and normalise, as in the non-conjugate slide. Any model, up to about three parameters, then the mesh explodes
Laplace and variational — replace the posterior by a normal (or another simple family) fitted near the mode. Fast, approximate
MCMC — simulate a Markov chain whose stationary distribution is the posterior; SSVS is one. Everything else, at the cost of convergence diagnostics

What changed in practice

None of this was usable for real models until cheap computation arrived. That is why the Bayesian–frequentist argument was philosophical for a century and is now largely a question of which tool fits the problem.

Where Economists Meet This

Macro and finance

DSGE estimation — short samples, many parameters, priors doing real work
Bayesian VARs — shrinkage is what makes them forecast
Model averaging — weight models by posterior probability instead of picking one

Micro and policy

Hierarchical models for regional and firm panels
Small-area estimation — official statistics for thin geographies
Sequential experiments — yesterday’s posterior is today’s prior

The recurring theme

Every one has more parameters than the data can comfortably identify. A prior is not a philosophical stance there — it is what makes the problem tractable.

Summary

Key Takeaways

The posterior is the whole answer; MAP, the posterior mean and the credible interval are all readings of it
Conjugate priors update by adding to their parameters: successes to \(\alpha\), failures to \(\beta\); precisions add in the normal case. Without conjugacy, compute the posterior numerically
The posterior mean shrinks. A normal prior shrinks everything alike; the horseshoe shrinks noise and spares large effects
Spike-and-slab priors turn variable selection into posterior inclusion probabilities
A credible interval licenses the probability statement people wrongly make about confidence intervals
A loss function is what turns a posterior into a decision — and defines “best”

Bayesian and frequentist

Agree asymptotically, under regularity conditions
Differ where it matters: small \(n\), many parameters, constrained spaces, and genuine prior information
Bayes factors can support \(H_0\); p-values never can

The cost, stated plainly

Every Bayesian answer is conditional on a prior. Report it, and show what happens when it changes.