{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "de2eb7e6",
   "metadata": {},
   "source": [
    "# Байесовский выбор моделей и информационные критерии"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31eb36fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import warnings\n",
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as mpatches\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.stats as st\n",
    "from scipy.special import comb, betaln\n",
    "from sklearn import linear_model\n",
    "from collections import Counter\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "sns.set_style(\"whitegrid\")\n",
    "sns.set_palette(\"colorblind\")\n",
    "palette = sns.color_palette()\n",
    "sns_colors = palette\n",
    "figsize = (12, 7)\n",
    "legend_fontsize = 16\n",
    "\n",
    "from matplotlib import rc\n",
    "rc('font', **{'family': 'sans-serif'})\n",
    "rc('text', usetex=True)\n",
    "rc('text.latex', preamble=r'\\usepackage[utf8]{inputenc}\\usepackage[russian]{babel}')\n",
    "rc('axes', **{'titlesize': '16', 'labelsize': '16'})\n",
    "rc('legend', **{'fontsize': '16'})\n",
    "rc('figure', **{'dpi': 200})\n",
    "\n",
    "SEED = 2026\n",
    "np.random.seed(SEED)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaf12fa0",
   "metadata": {},
   "source": [
    "# 1. Обоснованность и бритва Оккама"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9dfa529f",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 20\n",
    "ks = np.arange(n + 1)\n",
    "p_k_M0 = comb(n, ks) * 0.5 ** n          # биномиальное распределение\n",
    "p_k_M1 = np.full(n + 1, 1 / (n + 1))     # при равномерном априорном все k равновероятны\n",
    "\n",
    "fig, ax = plt.subplots(figsize=figsize)\n",
    "ax.bar(ks - 0.2, p_k_M0, width=0.4, label=r'$\\mathcal{M}_0$: честная монетка', color=palette[0])\n",
    "ax.bar(ks + 0.2, p_k_M1, width=0.4, label=r'$\\mathcal{M}_1$: $\\theta\\sim\\mathrm{Beta}(1,1)$', color=palette[1])\n",
    "ax.set_xlabel('число орлов $k$ из $n=20$')\n",
    "ax.set_ylabel(r'$p(k\\mid\\mathcal{M})$')\n",
    "ax.legend(fontsize=legend_fontsize)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b4dbccc",
   "metadata": {},
   "outputs": [],
   "source": [
    "log_ev_M0 = n * np.log(0.5) * np.ones(n + 1)\n",
    "log_ev_M1 = betaln(ks + 1, n - ks + 1)\n",
    "log_bf = log_ev_M1 - log_ev_M0\n",
    "\n",
    "fig, ax = plt.subplots(figsize=figsize)\n",
    "ax.plot(ks, log_bf, marker='o', linewidth=2, color=palette[2])\n",
    "ax.axhline(0, color='black', linewidth=1)\n",
    "ax.set_xlabel('число орлов $k$ из $n=20$')\n",
    "ax.set_ylabel(r'$\\log\\frac{p(D\\mid\\mathcal{M}_1)}{p(D\\mid\\mathcal{M}_0)}$')\n",
    "plt.show()\n",
    "\n",
    "for k in [10, 14, 15, 16, 18]:\n",
    "    print(\"k = %2d: байесовский фактор в пользу M1 = %8.3f\" % (k, np.exp(log_bf[k])))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f546cdd",
   "metadata": {},
   "source": [
    "# 2. Обоснованность в полиномиальной регрессии"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32f0d8f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "orig = lambda x: np.sin(2 * x)\n",
    "xd = np.arange(-1.5, 2, 7 / 60)\n",
    "num_points = len(xd)\n",
    "data = orig(xd) + np.random.normal(0, .25, num_points)\n",
    "xs = np.arange(xd[0] - 1.5, xd[-1] + 1.5, 0.01)\n",
    "\n",
    "xs_d = np.vstack([xs ** i for i in range(1, num_points + 1)]).transpose()\n",
    "xd_d = np.vstack([xd ** i for i in range(1, num_points + 1)]).transpose()\n",
    "\n",
    "\n",
    "def get_X(d, xx=xd_d):\n",
    "    return np.hstack([np.ones(xx.shape[0]).reshape(-1, 1), xx[:, :d]])\n",
    "\n",
    "\n",
    "def bayesian_update(mu, Sigma, X, t, beta=4):\n",
    "    Sigma_inv = np.linalg.inv(Sigma)\n",
    "    Sigma_n_inv = Sigma_inv + beta * X.T.dot(X)\n",
    "    Sigma_n = np.linalg.inv(Sigma_n_inv)\n",
    "    mu_n = Sigma_n.dot(Sigma_inv.dot(mu.T) + beta * X.T.dot(t))\n",
    "    return mu_n, Sigma_n\n",
    "\n",
    "\n",
    "def log_marginal_likelihood(X, t, alpha, beta):\n",
    "    N, M = X.shape\n",
    "    mu, sigma = bayesian_update(np.zeros(M), (1 / alpha) * np.identity(M), X, t, beta=beta)\n",
    "    A = np.linalg.inv(sigma)\n",
    "    return 0.5 * (M * np.log(alpha) + N * np.log(beta) - np.linalg.slogdet(A)[1]\n",
    "                  - beta * np.sum((t - X.dot(mu)) ** 2) - alpha * np.sum(mu ** 2) - N * np.log(2 * np.pi))\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=figsize)\n",
    "ax.set_xlim((xs[0], xs[-1])); ax.set_ylim((-2, 2))\n",
    "ax.scatter(xd, data, marker='*', s=120, color='black', zorder=5)\n",
    "ax.plot(xs, orig(xs), linewidth=1, label='Исходная функция', color='black')\n",
    "for d in [1, 3, 7]:\n",
    "    X = get_X(d)\n",
    "    mu, _ = bayesian_update(np.zeros(d + 1), (1 / 0.01) * np.identity(d + 1), X, data, beta=4)\n",
    "    ax.plot(xs, get_X(d, xs_d).dot(mu), linewidth=2, label='$d=%d$' % d)\n",
    "ax.legend(loc='lower center', fontsize=legend_fontsize)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72692372",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = np.arange(1, 12)\n",
    "fig, ax = plt.subplots(figsize=(12, 5))\n",
    "for alpha in [0.005, 0.01, 0.1, 1]:\n",
    "    ax.plot(ds, [log_marginal_likelihood(get_X(d), data, alpha, 4) for d in ds], marker='o',\n",
    "            label=r'$\\alpha=%g$' % alpha)\n",
    "ax.set_xlabel('степень многочлена $d$')\n",
    "ax.set_ylabel(r'$\\ln p(D\\mid d, \\alpha, \\beta=4)$')\n",
    "ax.set_xlim((1, 11))\n",
    "ax.legend(fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47b97120",
   "metadata": {},
   "source": [
    "# 3. Эмпирический байес"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4b1a533",
   "metadata": {},
   "outputs": [],
   "source": [
    "alphas_grid = np.logspace(-5, 5, 100)\n",
    "fig, ax = plt.subplots(figsize=(12, 5))\n",
    "ax.set_xscale('log')\n",
    "for i, beta in enumerate([1, 4, 10]):\n",
    "    ax.plot(alphas_grid, [log_marginal_likelihood(get_X(3), data, a, beta) for a in alphas_grid],\n",
    "            color=sns_colors[i], label=r'$d=3$, $\\beta=%g$' % beta)\n",
    "    ax.plot(alphas_grid, [log_marginal_likelihood(get_X(10), data, a, beta) for a in alphas_grid],\n",
    "            color=sns_colors[i], linestyle='--', label=r'$d=10$, $\\beta=%g$' % beta)\n",
    "ax.set_xlim((1e-5, 1e5)); ax.set_ylim((-50, -10))\n",
    "ax.set_xlabel(r'$\\alpha$')\n",
    "ax.set_ylabel(r'$\\ln p(D\\mid\\alpha,\\beta)$')\n",
    "ax.legend(fontsize=14, ncol=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "383e86a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def empirical_iteration(X, t, XtX, XtT, lambdas, alpha, beta):\n",
    "    N, d = X.shape\n",
    "    cur_lambdas = beta * lambdas\n",
    "    gamma = np.sum(cur_lambdas / (alpha + cur_lambdas))\n",
    "    A = beta * XtX + alpha * np.identity(d)\n",
    "    mu = beta * np.linalg.inv(A).dot(XtT)\n",
    "    return gamma / np.sum(mu ** 2), (N - gamma) / np.sum((t - X.dot(mu)) ** 2)\n",
    "\n",
    "\n",
    "def fit_eb(X, t, max_iter=400, eps=1e-4):\n",
    "    XtX, XtT = X.T.dot(X), X.T.dot(t)\n",
    "    lambdas = np.linalg.eigvalsh(XtX)\n",
    "    alpha, beta = [1.], [1.]\n",
    "    for _ in range(max_iter):\n",
    "        a, b = empirical_iteration(X, t, XtX, XtT, lambdas, alpha[-1], beta[-1])\n",
    "        alpha.append(a); beta.append(b)\n",
    "        if np.isclose(alpha[-1], alpha[-2], eps) and np.isclose(beta[-1], beta[-2], eps):\n",
    "            break\n",
    "    return alpha, beta\n",
    "\n",
    "\n",
    "alphas3, betas3 = fit_eb(get_X(3), data)\n",
    "alphas10, betas10 = fit_eb(get_X(10), data)\n",
    "print(\"d = 3:  alpha = %.4f, beta = %.2f (шум sigma = %.3f), итераций %d\" % (alphas3[-1], betas3[-1], betas3[-1] ** -.5, len(alphas3) - 1))\n",
    "print(\"d = 10: alpha = %.4f, beta = %.2f (шум sigma = %.3f), итераций %d\" % (alphas10[-1], betas10[-1], betas10[-1] ** -.5, len(alphas10) - 1))\n",
    "print(\"настоящий шум: sigma = 0.25, beta = 16\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddf5b307",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 5))\n",
    "ax.set_xscale('log')\n",
    "for i, (d, al, be) in enumerate([(3, alphas3, betas3), (10, alphas10, betas10)]):\n",
    "    ev = [log_marginal_likelihood(get_X(d), data, a, be[-1]) for a in alphas_grid]\n",
    "    ax.plot(alphas_grid, ev, color=sns_colors[i], label=r'$d=%d$, $\\beta=%.2f$' % (d, be[-1]))\n",
    "    ax.axvline(al[-1], color=sns_colors[i], linestyle='--')\n",
    "ax.set_xlim((1e-5, 1e5)); ax.set_ylim((-50, 0))\n",
    "ax.set_xlabel(r'$\\alpha$'); ax.set_ylabel(r'$\\ln p(D\\mid\\alpha,\\beta)$')\n",
    "ax.legend(fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82ab34dd",
   "metadata": {},
   "source": [
    "# 4. BIC, AIC и обоснованность по степени многочлена"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a228bf5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def linregr_ml(X, y):\n",
    "    mu = np.linalg.lstsq(X, y, rcond=None)[0]\n",
    "    return mu, np.mean((X.dot(mu) - y) ** 2)\n",
    "\n",
    "\n",
    "N = len(data)\n",
    "ds = np.arange(1, 12)\n",
    "minus2logL, AICs, BICs, minus2ev = [], [], [], []\n",
    "for d in ds:\n",
    "    X = get_X(d)\n",
    "    _, s2 = linregr_ml(X, data)\n",
    "    m2l = N * (np.log(2 * np.pi) + np.log(s2) + 1)\n",
    "    minus2logL.append(m2l)\n",
    "    AICs.append(m2l + 2 * (d + 2))\n",
    "    BICs.append(m2l + np.log(N) * (d + 2))\n",
    "    al, be = fit_eb(X, data)\n",
    "    minus2ev.append(-2 * log_marginal_likelihood(X, data, al[-1], be[-1]))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=figsize)\n",
    "ax.plot(ds, minus2logL, marker='o', linewidth=2, label=r'$-2\\ell(\\mathbf{w}_{\\mathrm{ML}})$')\n",
    "ax.plot(ds, AICs, marker='o', linewidth=2, label='AIC')\n",
    "ax.plot(ds, BICs, marker='o', linewidth=2, label='BIC')\n",
    "ax.plot(ds, minus2ev, marker='o', linewidth=2, label=r'$-2\\ln p(D\\mid d)$, эмпирический байес')\n",
    "ax.set_xlabel('степень многочлена $d$')\n",
    "ax.set_xlim((1, 11))\n",
    "ax.legend(fontsize=legend_fontsize)\n",
    "plt.show()\n",
    "\n",
    "for name, vals in [('AIC', AICs), ('BIC', BICs), ('-2 ln p(D|d)', minus2ev)]:\n",
    "    print(\"%-14s выбирает d = %d\" % (name, ds[int(np.argmin(vals))]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1447d62",
   "metadata": {},
   "source": [
    "# 5. Откуда в AIC берётся $2d$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac3b4188",
   "metadata": {},
   "outputs": [],
   "source": [
    "true_w, sigma = -.5, .5\n",
    "orig_lin = lambda x: true_w * x\n",
    "xd5 = np.array([-3, -1, 0, 0.5, 2.5]) / 2\n",
    "data5 = np.array([0.66, 0.38, -0.27, -0.21, 0.23])\n",
    "\n",
    "\n",
    "def log_likelihood(w, data=data5):\n",
    "    return - np.sum((data - w * xd5) ** 2) / (2 * sigma ** 2)\n",
    "\n",
    "\n",
    "def expected_log_l(w):\n",
    "    return - len(xd5) / 2. - np.sum(((w - true_w) * xd5) ** 2) / (2 * sigma ** 2)\n",
    "\n",
    "\n",
    "def w_ml(data):\n",
    "    return np.sum(xd5 * data) / np.sum(xd5 ** 2)\n",
    "\n",
    "\n",
    "def plot_two_curves(data, xlim=(-1, 0.55)):\n",
    "    w_hat = w_ml(data)\n",
    "    ww = np.arange(xlim[0], xlim[1], 0.01)\n",
    "    fig, ax = plt.subplots(figsize=(10, 6))\n",
    "    ax.plot(ww, [log_likelihood(w, data) for w in ww], linewidth=2, color=palette[0],\n",
    "            label=r'$\\log p(\\mathbf{y}\\mid w)$ на выборке')\n",
    "    ax.plot(ww, [expected_log_l(w) for w in ww], linewidth=2, color=palette[1],\n",
    "            label=r'$\\mathrm{E}_{\\mathbf{y}\\sim p_{\\mathrm{data}}}\\left[\\log p(\\mathbf{y}\\mid w)\\right]$')\n",
    "    for w in (true_w, w_hat):\n",
    "        ax.axvline(w, color='black', linestyle='--', linewidth=1)\n",
    "    top = ax.get_ylim()[1]\n",
    "    ax.text(true_w, top, r'$w_{\\mathrm{true}}$', ha='center', va='bottom', fontsize=14)\n",
    "    ax.text(w_hat, top, r'$w_{\\mathrm{ML}}$', ha='center', va='bottom', fontsize=14)\n",
    "    ax.set_xlabel('$w$')\n",
    "    ax.legend(fontsize=14, loc='lower center')\n",
    "    plt.show()\n",
    "    print(\"w_ML = %.3f; на выборке log p = %.3f, в ожидании по новым данным = %.3f\"\n",
    "          % (w_hat, log_likelihood(w_hat, data), expected_log_l(w_hat)))\n",
    "\n",
    "\n",
    "plot_two_curves(data5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4920703c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_B(data):\n",
    "    w_hat = w_ml(data)\n",
    "    l_hat, l_0 = log_likelihood(w_hat, data), log_likelihood(true_w, data)\n",
    "    e_hat, e_0 = expected_log_l(w_hat), expected_log_l(true_w)\n",
    "    return [l_hat - l_0, l_0 - e_0, e_0 - e_hat]\n",
    "\n",
    "\n",
    "B = np.array([get_B(orig_lin(xd5) + np.random.normal(0, sigma, len(xd5))) for _ in range(10000)])\n",
    "for i in range(3):\n",
    "    print(\"B%d = %.3f +- %.3f\" % (i + 1, B[:, i].mean(), B[:, i].std()))\n",
    "print(\"b = B1 + B2 + B3 = %.3f  (теория: d = 1)\" % B.sum(axis=1).mean())\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
    "for i, ax in enumerate(axes):\n",
    "    sns.histplot(B[:, i], ax=ax, kde=True, stat='probability', color=palette[i])\n",
    "    ax.set_title('$B_%d$' % (i + 1)); ax.set_ylabel(None)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}