{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3.4 贝叶斯模型的比较"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们从贝叶斯的角度来探讨模型选择的问题。\n",
    "\n",
    "假设我们希望比较 $L$ 个模型 $\\{\\mathcal M_i\\}, i=1,\\dots,L$。\n",
    "\n",
    "这里不同的模型表示在同一个观测集 $\\mathcal D$ 下的不同概率分布。\n",
    "\n",
    "假设第 $i$ 个模型的先验分布为 $p(\\mathcal M_i)$，则其后验分布为：\n",
    "\n",
    "$$\n",
    "p(\\mathcal M_i|\\mathcal D) \\propto p(\\mathcal M_i)p(\\mathcal D|\\mathcal M_i)\n",
    "$$\n",
    "\n",
    "如果这 $L$ 个模型的先验分布是一样的，那么对于后验分布，我们只需要关心\n",
    "`model evidence` 这个量： $p(\\mathcal D|\\mathcal M_i)$（它也可以看成模型的边际似然函数，因为它可以看成是模型空间的一个似然）。通常，我们把不同模型的 `model evidence` 的比值 $p(\\mathcal D|\\mathcal M_i)/p(\\mathcal D|\\mathcal M_j)$ 叫做 `Bayes factor`。\n",
    "\n",
    "当我们知道不同模型的后验分布后，预测值的分布可以写成它们的混合（模型平均）：\n",
    "\n",
    "$$\n",
    "p(t|\\mathbf x,\\mathcal D) = \\sum_{i=1}^L p(t|\\mathbf x,\\mathcal M_i, \\mathcal D) p(\\mathcal M_i|\\mathcal D)\n",
    "$$\n",
    "\n",
    "在这种情况下，如果两个模型分别是 $t=a, t=b$ 处的单峰模型，那么它们的混合是一个 $t=a, t=b$ 的双峰模型。\n",
    "\n",
    "对这个混合分布的一个简单近似是使用单一的最好模型，这种方法叫做模型选择（`model select`）。\n",
    "\n",
    "给定模型的一组参数 $\\bf w$，其 `model evidence` 为：\n",
    "\n",
    "$$\n",
    "p(\\mathcal D|\\mathcal M_i) = \\int p(\\mathcal D|\\mathbf w, \\mathcal M_i) p(\\mathbf w|\\mathcal M_i) d\\mathbf w\n",
    "$$\n",
    "\n",
    "Bayes 理论给出\n",
    "\n",
    "$$\n",
    "p(\\mathbf w | \\mathcal D, \\mathcal M_i) = \\frac{p(\\mathcal D|\\mathbf w, \\mathcal M_i) p(\\mathbf w|\\mathcal M_i)}{p(\\mathcal D|\\mathcal M_i)}\n",
    "$$\n",
    "\n",
    "为了简单，我们考虑单参数 $w$ 的模型，参数的后验分布正比于 $p(\\mathcal D|w)p(w)$，这里我们忽略了 $\\mathcal M_i$，让表示更加简洁。\n",
    "\n",
    "假设后验分布在 $w_{MAP}$ 附近是一个锐峰，宽度为 $\\Delta_{posterior}$，我们可以将后验分布的积分近似为峰值乘以宽度，更进一步，我们假设先验分布是均匀的，即 $p(w)=1/\\Delta w_{prior}$，那么我们可以这样近似 `model evidence` $p(\\mathcal D)$ 如下：\n",
    "\n",
    "$$\n",
    "p(\\mathcal D)=\\int p(\\mathcal D|w)p(w)dw \\simeq p(D|w_{WAP}) \\frac{\\Delta w_{posterior}}{\\Delta w_{prior}}\n",
    "$$\n",
    "\n",
    "取对数 \n",
    "\n",
    "$$\n",
    "\\ln p(\\mathcal D) \\simeq \\ln p(D|w_{WAP}) + \\ln \\left( \\frac{\\Delta w_{posterior}}{\\Delta w_{prior}} \\right)\n",
    "$$\n",
    "\n",
    "近似的示意图如下。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADyCAYAAADutRY4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHshJREFUeJzt3Xt8TXe+//FXrtv9UgkRcaSOaDVMK6iWojo6NDqO86Da\nFFNt9WrwO50o0xs1KKXKmJ4z6BwqrSFz0AraQydxiWuNTgxpxwwal0xJ8piEJhq5rN8f6ySkSdh7\nZ2etvbPfz8djPzz2zf7sXN77k+/6ru83wDAMRETEGoF2FyAi4k8UuiIiFlLoiohYSKErImIhha6I\niIUUuiIiFlLoiohYKPgm92sSr3jUihUrWLt2LTt37rS7FJH6FFDrHTc5OUKhKx5TXFxMbGwsp0+f\nZtu2bQwdOtTukkTqi0JX7Ld06VK6du3KuHHj6NKlCwcPHrS7JJH6UmvoakxXLHHlyhUyMjJ46KGH\nmDZtGocPHyYlJcXuskQsp05XLPHOO+8wePBg4uLiKCoqonPnzkRGRnLkyJFqj/3kk084duwY+/bt\nY9OmTYSGhgIQFxfHmjVr6N69O4mJiZw7d45169ZZ/VZEnKFOV+xTWFjIiRMniIuLA6BJkybMmDGD\njIwMNmzYUOWxFy5c4NSpU8yYMYM9e/aQmZkJwPHjx8nIyCAsLAyA6OhoMjIyrH0jIh6g0JV6t2zZ\nMqZMmVLlthdeeIHIyEhmzZpV5fbt27eTkJBAWloaDoeD7t27A5Cenk5MTAwREREATJo0iW7dullS\nv7Nee+01Nm/ebHcZ4uUUulKvLl26xLlz54iNja1yu8Ph4NVXXyUzM7PKEMH48eOJiIhg7dq1JCQk\nEBxszmpMT09n0KBBlY8rLy+nX79+Hq93yJAhXLlyxa3nzpkzhxEjRni4ImloFLpSr5YtW8bUqVNr\nvG/ixIl06tSJN998kx8eW9iyZQujRo2qvL53714GDhxYef3AgQMMGDDAo7V+++23XL16lcaNG3v0\n/xW5nkJX6k1+fj55eXnExMTUeH9wcDBvvPEGJ06cICkpqfL2nJwccnNz6dmzJwAlJSVkZWVVGU5I\nT0+nb9++HDlyhBUrVjB9+nQ2btzI8uXLWbNmDQDZ2dnMnj2bzZs388Ybb3Dq1CnA7L4XL15MSkoK\nixcvBiA1NZVf/OIXdOzYkQ8//BAwx5FffvllUlJSmD17NgAHDx7kgw8+4Mc//jGLFy+mZ8+enD17\nlrVr1zJmzJjK+mp67Zqee+7cOU99ucVH3OyMNBG3LVy4EMMwWLp0aa2PKSsrIzQ0lNmzZzN27FiC\ngoJo3bo1TZs2pbCwkBYtWrBr1y4cDgelpaUAnDx5kg4dOgBw8eJFYmNjSU1NZcGCBRQVFXHnnXcy\nevRoRowYwWeffUZYWBghISHMnz+fFStW8NFHH5Gbm8vTTz9dOQb7wAMPsGrVKn7+85/Tt29fcnJy\niI+P59ChQ7Rr144DBw4A5kHAbt26ERoayksvvcSkSZPYs2cP8fHxLFmyBICioqIaX3vy5MnVnutw\nOOrzWyBeSKEr9aKgoIBly5ZRWFjo1ONPnz5NUlISEyZMIDg4mKSkJBITE7n99tsJDw8nOTmZ+fPn\n07t3b5o3b87kyZMBGDZsGDNnzuThhx8G4Msvv6R169asX7+euLi4ytkOX331VeWwwciRI/nlL39J\n9+7def311ytrOHToEH369AHgD3/4A9HR0Rw9epSLFy/y4osvAtCjRw/mzJlT2dU6HA6GDBnCu+++\nyxNPPAFQ62vX9FzxQ4Zh3Ogi4vXuuece48yZM4ZhGMbEiRON5ORkY/ny5cbbb79tGIZhXLlyxbj3\n3nuNM2fOGKmpqcb06dMNwzCMf/zjH8Zjjz1mGIZhHD9+3PjpT39qGIZhrFu3zli9erUxd+7cytc4\nceKEUVRUZBiGYQwcONA4f/58lRri4uKMgoICIyUlpcbXPnv2bK3PlQap1lxVpys+LT8/n/z8fNLS\n0iguLubuu+/mkUce4dKlS7z11lts3bqVo0eP8v7779OxY0cuX75MTEwMmzdv5ptvvmHRokUAtGnT\nhltuuYX169fTv39/wsLC+NWvfsWWLVu4evUqLVu2JCYmhvLycgIDA4mMjKxSR0xMDCkpKQwbNoyQ\nkJBqrx0VFVXrc8W/6Iw08WkbNmzg0KFDLFiwwO5SRK6nM9Kk4fn6669ZsmQJeXl5fPfdd3aXI+IU\ndboiIp6nTldExBsodMUj8vPzmT17doMbW920aRMvvvgiWVlZdpciDYRCV+qkImy7dOnCqVOnqpyV\n1RAMHjyYVq1aERcXx/PPP6/wlTrTmK647dChQ/Tr1w+Hw0F0dDSNGjWyu6R6U1paypkzZygoKGD+\n/Pm8/PLLdpck3q3WMV3N0xW39erVi9dee42VK1fSrFkznn32We666y67y/K4rKwsfve733HmzBme\ne+45Jk6caHdJ4sPU6UqdlZSUsGbNGubMmUNCQgLz5s2zuySPSU5OZtKkSUyZMoUpU6bQsmVLu0sS\n36CNKaX+lZSUcOnSJdq0aWN3KR5TWFhIeXk5zZs3t7sU8S0KXRERC2meroiIN1DoiohYSKErImIh\nha6IiIUUuiIiFlLoiohYSKErImIhha6IiIW09oLYwzDg0iUoL7fn9Vu2hED1HGI9ha5YwzBg3z74\n9FM4csS8fPcdhIRYX0tF0PfsCXFxcP/9MHw4BAVZX4v4HZ0GLPWrrAw+/hgWLYKcHEhIgN69zbCL\nioKAWs+WrF+5ufDll2b4b9pkXn/pJZgwAZo0sacmaUi09oLY4LvvoH9/aNoUEhPh3/7NO7vJii58\n0SLYvx/27IGYGLurEt+m0BUbPP88FBfDqlV2V+K8Zctg7VozeIM1+iZu04I3YrFPP4XPPoOlS+2u\nxDWTJkGzZvD223ZXIg2UOl3xvLw8+NGP4KOPzINUvubsWejVC/73f82DbSKuU6crFjEMePFFePRR\n3wxcgI4dYfFiGD8evv/e7mqkgVHoimft3m3OCvD1LXvGjjUPpv3nf9pdiTQwCl3xrM8/h0ceAV/f\nGTggAMaNgz/+0e5KpIFR6Ipn7dkDAwbYXYVn3HefOZWsrMzuSqQBUeiK5xQXw+HD0K+f3ZV4Rrt2\n0LYtHDtmdyXSgCh0xXP+9Ce47TZo0cLuSjxn4ECzexfxEIWueE5DGlqoMGCAQlc8SqErnrN7d8MM\n3d27zalwIh6g0BXPKCszDzrdd5/dlXhWdLS5XsTJk3ZXIg2EQlc849gx86BTu3Z2V+JZAQEaYhCP\nUuiKZzTE8dwKCl3xIIWueIZCV8QpCl2pO8No2KEbG2su4vPtt3ZXIg2AQlfq7tQpc+zz1lvtrqR+\nBAaai7Gr2xUPUOhK3VV0uXZtvWMFDTGIhyh0pe4a4lSxHxowAPbutbsKaQAUulJ3p06Zp/82ZF27\nmu9TpI4UulJ3586ZO/s2ZLfcYi5o/t13dlciPk6hK3VjGP4RugEB5ns8f97uSsTHKXSlbgoKzNNk\nmze3u5L6FxVlfsCI1IFCV+rGH7rcCgpd8QCFrtSNQlfEJQpdqZvz56FDB7ursEaHDhrTlTpT6Erd\nqNMVcYlCV+pGoSviEoWu1I1CV8QlCl2pG38K3bZtIT/fPElCxE0KXakbfwrdwECIjITsbLsrER+m\n0BX3FRZCcTG0bm13JdbREIPUkUJX3Hf+vBlCDXlJxx/SqcBSRwpdcd+5c/4zR7dChw7qdKVOFLri\nPn8az62g4QWpI4WuuE+hK+Iyha64T6Er4jKFrrhPoSviMoWuuM8fQzciAnJyoKTE7krERyl0xX3+\nGLohIRAeDt9+a3cl4qMUuuKe4mJz14jwcLsrsZ6WeJQ6UOiKe7KzoX1789RYf6NxXakDP/yNEY/w\nx6GFCgpdqQOFrrhHoWt3FeKjFLriHoWu3VWIj1LoinsUunZXIT5KoSvuUejaXYX4KIWuuKdiWUd/\nFBkJ//gHlJfbXYn4IIWuuMcfl3Ws0KgRtGhhnpkm4iKFrriutBQuXjRPifVXUVFw9qzdVYgPUuiK\n6/LyzC16QkLsrsQ+ERHmB4+IixS64rqLF/3z9N/rhYcrdMUtCl1xXU6OuR25P2vbVmO64haFrrhO\nna46XXGbQldcp05Xna64TaErrlOnq05X3KbQFdep01WnK25T6Irr1Omq0xW3KXTFdep0r3W6hmF3\nJeJjFLriOnW60LSp+W9hob11iM9R6Irr1OmaNK4rblDoimtKSuDyZWjVyu5K7KdxXXGDQldck5sL\nbdr454aUP6ROV9yg3xxxzcWLGlqooE5X3KDQFdfk5OggWgV1uuIGha64Rp3uNep0xQ0KXXGNOt1r\n1OmKGxS64hp1uteo0xU3KHTFNTox4hp1uuIGha64RidGXKNOV9yg0BXXqNO9Jjxc6y+IyxS64hp1\nutc0aQLBweYZeiJOUuiKa9TpVqVxXXGRQlecV1wMV65o3YXrVQwxiDhJoSvOy8mBsDAICLC7Eu/R\ntq0OpolLFLriPI3nVqdOV1yk0BXnaTy3OnW64iKFrjhPnW516nTFRQpdcZ463erU6YqLFLriPHW6\n1anTFRcpdMV56nSrU6crLlLoivPU6VanTldcpNAV56nTrU7rL4iLFLriPHW61TVqZF4uXbK7EvER\nCl1xnjrdmmmJR3GBQlecc+UKXL0KLVrYXYn30aI34gKFrjinYm80rbtQnTpdcYFCV5yj8dzaqdMV\nFyh0xTkaz62dOl1xgUJXnKNOt3Y6QUJcoNAV55w/D+3b212Fd2rfHrKz7a5CfIRCV5yTlQXR0XZX\n4Z2io82vj4gTFLrinG++gU6d7K7CO3XqZH59RJyg0BXnqNOtXbt25o7ARUV2VyI+QKErN2cYZuiq\n061ZQAD8y79oiEGcotCVm8vJgSZNoFkzuyvxXhrXFScpdOXmNJ57cxrXFScpdOXmNJ57c+p0xUkK\nXbk5dbo3p05XnBRsdwFWKi+H0lIoKzOPDVVcQGtQ34jjb1mUd+lKyWW7K/FeQW2jcZzKokhfo1pV\nrJUUEHDtEhQEwcEQ6EftX4MJ3fx82LEDvvoKMjPhb38zb7t82bwUF5vBGhxsfqOv/8Zf/8Mg1f3h\nyjesDnmQrb+yuxLvFVneiZ1XvqFrpN2VeKfrm5vrL+XlUFJi3hcaCs2bm8drW7SAzp2hWze44w64\n/37o2NG28j0qwLhxi+f1/d+lS7B0qXm55x740Y/Mb9Jtt8Ett5jfwGbNoHFj//o09agePeDDD+HO\nO+2uxHuVlUHTplBQAA6H3dX4nPJyszEqLDSbpIIC+PvfzQYqM9NsqMaMgVde8ZnwrbWF8+lOd+VK\nePVVGDYM9u+HmBi7K2qADENjus4ICoIOHeDMGf0guiEw0GyMGjeGsDDztrvuunZ/Tg4sWmTe9uST\n8PbbvttE+WjZsHMnzJwJu3bBmjX6Oa83//ynGSitWtldiffTDIZ6Ex4OCxaYw4d79sCvf213Re7z\nyU43Lw/Gj4dVq8wxH6lH6nKdpxkM9a5tW/j9782hxEGDoGdPuytync91uoYBTz8Njz4KQ4faXY0f\n0Bxd56nTtUTnzrBkCSQkmGPAvsbnQve//gvOnoV58+yuxE+o03WeOl3LPP642e1OnWp3Ja7zqdA9\nfx7eeMP88yI01O5q/IQ6Xeep07XUsmXmMZ0dO+yuxDU+Fbpr1sCoUdC1q92V+BF1us5Tp2up5s0h\nMRFWrLC7Etf4TOgaBqxeDRMm2F2Jn1Gn67yoKLhw4dpsf6l3jz5qdrp5eXZX4jyfCd39+80zxu65\nx+5K/Iw6XecFB0NEBJw7Z3clfqNVK4iPN4ccfYXPhO6qVWaXq1N1LVRQYHZtbdrYXYnv0Liu5SZM\nMPPBV/jEPN2iIvif/4Fjx+yuxM9U7BahTzrnaVzXcj/+MVy8CEePmssAeDuf6HQ3bjSHFTp0sLsS\nP6PxXNep07VcUBD87GfmMR9f4BOhu2qVeb61WEzjua5Tp2uLCRPgo4984xim14duVhZkZMCIEXZX\n4ofU6bpOna4tYmLMqaTbttldyc15feiuWWMu6daokd2V+CF1uq5Tp2ubCRN8Y4jBq0PXMMw/GcaP\nt7sSP/XXv8K//qvdVfiWjh3h22/No79iqdGjITXVXBjPm3l16H75JVy9qrm5tjh/HrKzqy5qKjfn\ncEDfvpCWZnclfqdlS3jwQdiwwe5KbsyrQ3ftWnNhC81YssGnn8JPfmJO+BfXxMfD1q12V+GXHn/c\nzA1v5rWhW1ZmnmXy+ON2V+Kntm0zw0NcN3y4+fXTbqeWi4+HP//Z/EPNW3lt6O7ebS5YfMcddlfi\nh4qL4Y9/NPdBEtdVrKyfmWlvHX6oUSP493+HdevsrqR2Xhu6FUMLYoP0dDM4wsPtrsQ3BQRc63bF\ncmPHevcQg1eGbnGxeRbaY4/ZXYmf2rbNDA1xn8Z1bTNokDmB5Ouv7a6kZl4Zup9+au767SNbLTc8\nW7dqPLeuBg+GI0fMRYPEUkFBZsPmrd2uV4auhhZsdPIk5Of75o5/3qRJE7jvPti+3e5K/FLFLAZv\nPJbpdaF79Ki5vfojj9hdiZ+qmLUQ6HU/Gr4nPl7jujaJi4PWrc2Tq7yNV/1mlZfD88/D3LnmF0xs\noKlinhMfb46VlZfbXYnfCQiA996DadO87ww1rwrd//5v89+nn7a3Dr9VWGjOXHjwQbsraRg6dza7\nh8OH7a7EL919tzl97JVX7K6kqgDjxoMelo2I5ORAbKy539Gdd1r1qlLJMMxFSQ0DPvzQ7moajnff\nNf/G3bULmja1uxq/k59vzvX/+GMzhC1U63m0XhO6Tz5pNgWLF1v1ilLFzJnw2WfmmgFNmthdTcNh\nGPDUU+bOiZs2mYfWxVIffmjmyqFDlp7VXmvo2j68UFoKCxfC55/Dm2/aXY2fWr0akpJg82afCtyd\nO3faXcLNBQTA8uXm0M1LL9ldjV8aO9bcwPK55+DSJbursTl0v/gC+vQxZ9Xs3GnuYy8WKi42w3b6\ndHNubrt2dlfkEp8IXYDQUHPpq88/N/+iyM+3uyK/EhBgfvkDA82hhk2b7K3H0iWkCgvN5RoPH4Z9\n+8z1FRYtMj+JtJKYRXJzza041q83fxJ79DB/CivWC5D60aqVOTPkP/7DXOj8wQfNyaR9+0JkpH4B\n6lnr1rBypZk5zz0Hv/mNuaFl797Qq5e1G17fcEz3i/D4Gu+s7caK241yKDfMf8vKobTE3LvIwOxm\nW7Y0L+3aQWhInd+DOKOgwFyUvLTU/LgfORISEnz6tL9Zs2Yxa9Ysu8tw3T//aX7gJSebE9MLC829\nZtq2VfhaoKwcLl6A/ALz1+JSgZlXIcEQEmIOuwcEmt+KwAAg4P8GaGv41tT23eqTs829A2kBAQFe\neD6HiIj3MwyjxuC94fDCTWY2iIhIzbx39oKIiD9R6IqIWEihKyJiIYWuiIiFFLoiIhZS6IqIWEih\nKyJiIYWuWOaTTz5h7ty5DB8+nKtXr1beHhcXx7FjxwBITEzkMe1IKg2YQlcsceHCBU6dOsWMGTPY\ns2cPmZmZABw/fpyMjAzCwsIAiI6OJiMjw85SReqVQlcssX37dhISEkhLS8PhcNC9e3cA0tPTiYmJ\nISIiAoBJkybRrVs3srKyGD58OM2aNWPt/23rumrVKhwOB1OnTuX8+fMAzJs3j169erH9BxtAlpSU\ncHcNq1bv27ePIUOGMHjwYFasWMF7773H888/z+7du+vz7YtcYxjGjS4iHvXkk08akydPrrw+btw4\n49lnn628XlpaaixcuNAwDMNYvXq1ceutt1bed/nyZaNJkybGn//858rbtmzZYuTn51d7nV//+tdG\nUFCQkZeXV+2++Ph4Izk5ufL6X/7yF6Ndu3Z1e2MiVdWaq+p0xVJbtmxh1KhRldf37t3LwIEDK68f\nOHCAAQMGAFQOOVRISkoiOjqa3NxcAPLy8ggODqZly5ZVHpeZmUlYWBi33HIL2dnZVe4rKysjPT2d\n/v37V952+vRpmvjQ4u3i2xS6YpmcnBxyc3Pp2bMnYA4BZGVl0e26tXzT09Pp27cvUDV09+3bR48e\nPWjfvn1l6G7bto2hQ4dWeQ3DMEhKSiIhIYF27dpVC90vvviCNm3aEBkZCcD333/PypUree+99zz/\nhkVqYOki5uLfWrduTdOmTSksLKRFixbs2rULh8NBaWkpACdPnqRDhw6Vj68I3atXr7J7925mzJhB\neHg4ubm5HDp0iHvuuafaayxfvpwnnngCgIiIiGqhm5qaSocOHUhOTubq1atcvnyZpUuXcuutt9bX\n2xapQqErlgkODiYpKYnExERuv/12wsPDSU5OZv78+fTu3ZvmzZszefLkyseHhYVhGAa/+c1vmDhx\nIgDh4eFkZ2dz8uRJEhISqvz/Z86cITU1lcDAQHbt2kVhYWG10E1LS+Opp55izJgx9f+GRWqg0BVL\njRw5kpEjR1a57eGHH67xsS1btiQ/P59WrVpVdr1t27Zl06ZNvP7669Ue/9vf/paPPvqIkBBzO5Kz\nZ89WCd3i4mL27dvHBx984Km3I+IyjemKV3vggQd46qmnKq9HRkYye/ZsGjVqVHnbgQMHGDFiBCdO\nnKi87eDBg/zpT39iz5497Nixg8OHDzNjxgxCQkLYsWOHpe9B5Ho33K6HWrZDExGRG9LOESIi3kCh\nKyJiIYWuiIiFFLoiIhZS6IqIWEihKyJiIYWuiIiFFLoiIhby2dDdv38/jzzyCH//+9/tLkVELFJS\nUsLIkSPZuHEj5eXldpfjFp8L3f379zNs2DAee+wxHnzwQa0OJeJHQkJCeOaZZ5g7dy49e/b0yfD1\nqdOAf/KTn7Bjxw7at29PeHg4gYE+95khIh5gGAYFBQVkZWURFRXFmTNn7C7ph2o9DdinVhlbtGgR\ns2bNIjU1lT59+jBu3LhquwaISMNWXFzMxx9/zKpVq+jfvz+vvPKK3SW5xKc63QrffPMN8+bNY8OG\nDRw8eJAuXbrYXZKIWKC0tJSuXbtyxx13MHPmTPr06WN3SbWptdP1ydCtkJeXR6tWrQgKCrK7FBGx\nyMWLF2nbtq3dZdxMwwxdEREvpaUdRUS8gUJXLLVixQruv/9+u8sQsY2GF8QyxcXFxMbGcvr06Rq3\nTxdpQDSmK/ZbunQpXbt2Zdy4cXTp0oWDBw/aXZJIfdGYrtjrypUrZGRk8NBDDzFt2jQOHz5MSkqK\n3WWJWE6drljinXfeYfDgwcTFxVFUVETnzp2JjIzkyJEj1R77ySefcOzYMfbt28emTZsIDQ0FIC4u\njjVr1tC9e3cSExM5d+4c69ats/qtiDhDna7Yp7CwkBMnThAXFwdAkyZNmDFjBhkZGWzYsKHKYy9c\nuMCpU6eYMWMGe/bsITMzE4Djx4+TkZFBWFgYANHR0WRkZFj7RkQ8QKEr9W7ZsmVMmTKlym0vvPAC\nkZGRzJo1q8rt27dvJyEhgbS0NBwOB927dwcgPT2dmJgYIiIiAJg0aRLdunVzu6bXXnuNzZs3u/18\nEXcpdKVeXbp0iXPnzhEbG1vldofDwauvvkpmZmaVIYLx48cTERHB2rVrSUhIIDjYXB4kPT2dQYMG\nVT6uvLycfv36uV3XnDlzGDFihNvPF3GXQlfq1bJly5g6dWqN902cOJFOnTrx5ptv8sNjC1u2bGHU\nqFGV1/fu3cvAgQMrrx84cIABAwbUT9Ei9cinVhkT35Kfn09eXh4xMTE13h8cHMwbb7zB008/TVJS\nEj/72c8AyMnJITc3l549ewLmwtVZWVlVhhPS09OZPn06R44c4fDhw5w8eZK+ffuSk5ND48aNue22\n2/j6669Zs2YNw4cPJykpiZSUFBo1asT27dv5+OOPSU5OBiA7O5v333+fu+66i8OHDzNhwgRycnJq\nfH5UVFQ9f9WkoVPoSr1ZuHAhhmGwdOnSWh9TVlZGaGgos2fPZuzYsQQFBdG6dWuaNm1KYWEhLVq0\nYNeuXTgcDkpLSwE4efIkHTp0AMzFT2JjY0lNTWXBggUUFRVx5513snHjRrp160ZoaCgvvfQSkyZN\nwuFw8PnnnxMfH8+SJUsAKCoqYsSIEXz22WeEhYUREhLC/PnzmTx5co3PF6krha7Ui4KCApYtW0Zh\nYaFTjz99+jRJSUlMmDCB4OBgkpKSSExM5Pbbbyc8PJzk5GTmz59P7969ad68OZMnTwZg2LBhzJw5\nk4cffhiAL7/8ktatW9OjRw/mzJnDmDFjACoDc8iQIbz77rs88cQTAKxfv564uLjKWRFfffUVjRs3\nrvX5InWlebri8+69916Sk5Pp2LEjzzzzDEOHDmX06NEMGjSI3//+90RGRlZ5fK9evUhLS2P37t1k\nZ2dTUFDAtGnT+P7773nggQdITk4mKiqq1ueLOKFh7Bwh8kP5+fnk5+eTlpZGcXExd999N6NHj6a8\nvJzAwMAaAzMmJoaUlBSGDRtGSEgIb731Flu3buXo0aO8//77REVF3fD5InWhTld82oYNGzh06BAL\nFiywuxSR67m94I2I1woICLgdWAn8Ffh/hmF8Z3NJIjel0BURsZBOjhARsZBCV0TEQgpdERELKXRF\nRCyk0BURsZBCV0TEQgpdERELKXRFRCyk0BURsdD/B8jEYodZRWTaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x6f4c1d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "_, ax = plt.subplots()\n",
    "\n",
    "xx = np.linspace(-1, 1, 100)\n",
    "yy_1 = np.ones_like(xx) * 2\n",
    "\n",
    "yy_1[:20] = 1 + np.tanh((xx[:20] + 0.8) * 30)\n",
    "yy_1[-20:] = 1 - np.tanh((xx[-20:] - 0.8) * 30)\n",
    "\n",
    "yy_2 = np.zeros_like(xx)\n",
    "\n",
    "yy_2[35:45] = 4 * (1 + np.tanh((xx[35:45] + 0.2) * 30))\n",
    "yy_2[-45:-35] = 4 * (1 - np.tanh((xx[-45:-35] - 0.2) * 30))\n",
    "\n",
    "yy_2[45:-45] = 4 * 2\n",
    "\n",
    "ax.plot(xx, yy_1, 'b', xx, yy_2, 'r')\n",
    "ax.set_ylim([-3, 10])\n",
    "\n",
    "ax.spines['right'].set_color('none')\n",
    "ax.spines['top'].set_color('none')\n",
    "ax.spines['left'].set_color('none')\n",
    "ax.xaxis.set_ticks_position('bottom')\n",
    "ax.spines['bottom'].set_position(('data',0))\n",
    "\n",
    "ax.set_yticks([])\n",
    "ax.set_xticks([0])\n",
    "ax.set_xticklabels(['$w_{WAP}$'], fontsize=\"xx-large\")\n",
    "\n",
    "\n",
    "ax.annotate(\"\",\n",
    "            xy=(-0.8, -1.8), xycoords='data',\n",
    "            xytext=(0.8, -1.8), textcoords='data',\n",
    "            arrowprops=dict(arrowstyle=\"<->\",\n",
    "                            connectionstyle=\"arc3\"), \n",
    "            )\n",
    "\n",
    "ax.text(-0.15, -2.8, r\"$\\Delta w_{prior}$\", fontsize=\"xx-large\")\n",
    "\n",
    "ax.annotate(\"\",\n",
    "            xy=(-0.2, 8.8), xycoords='data',\n",
    "            xytext=(0.2, 8.8), textcoords='data',\n",
    "            arrowprops=dict(arrowstyle=\"<->\",\n",
    "                            connectionstyle=\"arc3\"), \n",
    "            )\n",
    "\n",
    "ax.text(-0.15, 9.3, r\"$\\Delta w_{posterior}$\", fontsize=\"xx-large\")\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于多参数变量，设参数个数为 $M$，并假设所有的参数有相同的 $\\frac{\\Delta w_{posterior}}{\\Delta w_{prior}}$，我们有\n",
    "\n",
    "$$\n",
    "\\ln p(\\mathcal D) \\simeq \\ln p(D|\\mathbf w_{WAP}) + M \\ln \\left( \\frac{\\Delta w_{posterior}}{\\Delta w_{prior}} \\right)\n",
    "$$\n",
    "\n",
    "当我们增加模型复杂度时，上式的第一项通常会减小，因为通常复制的模型对数据的拟合更好，但是第二项会随着 $M$ 的增大而增大；而我们的目标是最大化模型的 `model evidence`，因此要对这两项做一个权衡。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}