vi post
This commit is contained in:
ritchie46 committed
1 parent
8f581f8b4b
commit
70f439ab8e
1 file changed
+525
@@ -0,0 +1,525 @@
|
|||||||
|
{
|
||||||
|
"cells": [
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 8,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"import numpy as np\n",
|
||||||
|
"import matplotlib.pyplot as plt\n",
|
||||||
|
"from scipy import stats"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 11,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"lengths = np.array([195, 185])"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 12,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"# lets create a grid of our two parameters\n",
|
||||||
|
"mu = np.linspace(150, 250)\n",
|
||||||
|
"sigma = np.linspace(0, 15)[::-1]\n",
|
||||||
|
"\n",
|
||||||
|
"mm, ss = np.meshgrid(mu, sigma) # just broadcasted parameters"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 13,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"name": "stderr",
|
||||||
|
"output_type": "stream",
|
||||||
|
"text": [
|
||||||
|
"/opt/miniconda3/lib/python3.7/site-packages/scipy/stats/_distn_infrastructure.py:1662: RuntimeWarning: divide by zero encountered in true_divide\n",
|
||||||
|
" x = np.asarray((x - loc)/scale, dtype=dtyp)\n"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"text/plain": [
|
||||||
|
"Text(0, 0.5, '$\\\\sigma$')"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"execution_count": 13,
|
||||||
|
"metadata": {},
|
||||||
|
"output_type": "execute_result"
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAATgAAAEMCAYAAABQlKrRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFfxJREFUeJzt3X2MZXV9x/HPZ2Znn3dh2V0QZZfFJ4waCji1Wq1tpVqkRvyjf0C0oZVkU5P6lFqFkNakf6GYqkmbmI1SMBJIixhNqxVqpcRUsQuyCi4K9XF1KyzWKuzj7Hz7x73UYc737J7LnXvOvb95v5LJ3vnNuff8fjN3v/ec7+/JESEAKNFU1xUAgFEhwAEoFgEOQLEIcACKRYADUCwCHIBiEeAAFIsAB6BYBDgAxVrRdQWa2LJlc+zYvr3raiySzABJJ4UMcZwkxXxSlhw7nxw3fzx/zezY7DzZcXX1TDkpSsqmks/Zqelmx0mSs2OT8zStT2159vz0yflrYsnc8/X7DkTE1pMdNxEBbsf27dr95TtHfp502lrtVLYhgsz8XLVs7lh+mrmj1bJjR6q1OfREtezgL/LXPFw9Ni07Uj1PbT0z00ngWbm6WrZqTbVs7fpKkddUyyRJazZUj83Os2JlUjaTv+Z08l8jC7pOgm5algc91wVYnJDXnfqDJsdxiwqgWAQ4AMWaiFvU5WHIW5U0V5fl0CQdT26bs7LsdvR4cntdW6fkNn46ec355NZxfpB0QUOjuB3kFnOscQUHoFgEOADFIsABKBYBDkCx6GQ4mbokctOEd9MBo/UvMMCxi2SJeqlmUG82Xq9hx0Od7N2VdWYM0kGy1Go7CYb9u2EccAUHoFgjC3C2r7f9iO37k5+923bY3jKq8wPAKK/gbpB08eJC29skvUbSD0d4bgAYXQ4uIu6yvSP50YckvUfSZ0Z17qcrmxdYv61iw4nbTcemDpveGSRnleXW0oG+yaDeNIc2wADcbI5nVp90XnBLeblhMfh3bLSag7P9Bkk/jog9bZ4XwPLUWi+q7bWSrpH02obH75S0U5K2b9s2wpoBKFWbV3DPkXSOpD22vy/pLEn32n5GdnBE7IqI2YiY3bplc4vVBFCK1q7gIuKbkk5/8vt+kJuNiANt1WFstJWjqctZpePbGq5Rlx03SA5upmH+L1szbxCNxx8OsODlEGMaWfetG6McJnKzpK9IOtf2PttXjupcAJAZZS/q5Sf5+Y5RnRsAJGYyACgYAQ5AsZhsv5QaJ6YbPneQ52fqkv/ZJPx0sn3SSXEs2QRnEE07KYZdvRcQV3AACkaAA1AsAhyAYpGDe7rSifXDLII5pEEmpzedbH80ybelObQBJsE33dEr3UR7kLxcw13oB8l9Mlh34nAFB6BYBDgAxSLAASgWAQ5AsehkOIm6VSDqV/qtvMJwxw010Lcm+Z+u3pusHHI8KRtkoK+Tz8/jSSdF0xV9hzWKTgI6HsYaV3AAikWAA1AsAhyAYpGDW1JDrAI7SC6n6bF1eawsN5eu6Ntsld+6fGSav0xzfQ13u68zzOq7A+1s37Q65OXGBVdwAIpFgANQLAIcgGKRg+vCQBO5h8jnZBPWpZrJ7UlZMuYtjiU5tLpcX5qDazgOrjWD/H7JrU2aUe6qdb3tR2zfv6DsOtsP2v6G7U/bPnVU5weAUd6i3iDp4kVld0h6cUScJ+k7kq4e4fkBLHMjC3ARcZekny0quz0inrxH+ap6u9sDwEh02cnwFkmf7/D8AArXSSeD7WskzUm66QTH7JS0U5K2b9vWUs1GYNhBn00HrA4iS+onyf/IOgSyVX5rRFJ3N13Rd4DJ9unA2mHKTlSOidL6FZztKyS9XtKb4gRLckTEroiYjYjZrVs2t1dBAMVo9QrO9sWS3ivptyPiYJvnBrD8jHKYyM2SviLpXNv7bF8p6W8lbZB0h+37bH90VOcHgJFdwUXE5Unxx0d1vrZleZ/0fju7Cx8krzZMLql2sn1Snk2Czwb1DjvQN1swM5tYP8hk+/zkzepT+/RmxzKxfrwxVQtAsQhwAIpFgANQLCbbj9wIFl4cNu+TbjqTjHnL8m3ZOLi6Sf1Tyedn4wUvB9h0ZiQLXjY8D8YaV3AAikWAA1AsAhyAYhHgABSLToZxUTvpu2Fhtot83WDZdFethoN6jxxJXq+mQ2B6unpo0qHguk6KxoaZWD+ClZQxNriCA1AsAhyAYhHgABSLHNwSSifgtzXZfhBJHizLjaX5tmygb10ObkXy9soGFDfd2X5qkMUpR/B7w8ThCg5AsQhwAIpFgANQLHJwozZ0fqjpGK8BPqvSBS+bTbaPw83HwTnLwWVj67IcXPqCNW3Mfh9pvm64vByLW04eruAAFIsAB6BYBDgAxRrlrlrX237E9v0Lyk6zfYfth/r/bhrV+QFglJ0MN6i3TeAnFpRdJemLEXGt7av63793hHWYHHUJ7GxV3GE7Hpp2MiSDeucPVctqx/nOJB0S2Xma7qBVuyDBML8jdrYv2ciu4CLiLkk/W1R8qaQb+49vlPTGUZ0fANrOwZ0REfslqf/v6S2fH8AyMradDLZ32t5te/ejBx7rujoAJlDbA31/avvMiNhv+0xJj9QdGBG7JO2SpNkLLxhgi6VJMOTCi1leLllcsla2wGSSG4skB3f8cDJQdz7/88yvqr69ptIcXLYgQdLGQQb6jmKRAkyctq/gPivpiv7jKyR9puXzA1hGRjlM5GZJX5F0ru19tq+UdK2k19h+SNJr+t8DwEiM7BY1Ii6v+dFFozonACzEZPsRG24RzLryhmPe6l6z4Ti4+STflo6DO56PY5taNVMtG8lk+6b5ukFynyjB2PaiAsCwCHAAikWAA1AsAhyAYtHJ0IUs+V83lDldrTb5XJpKBvqmE/WVT25PBv/GkaTj4YnmnQzza5IduJpOts/aXTeYuWkHC6v3LjtcwQEoFgEOQLEIcACKRQ5uXAyymGNyrJN8WwySR0pyY8eTQb3HDlYXsayZa6/pbGL+MIN6szyjlOcamYAPcQUHoGAEOADFIsABKBYBDkCx6GToQL7CyACriTRdPWNqgD9vtprIkWonweHD1U6C4zW9DCsPNhzom0kHLtcN9B1gRzEsK7wLABSLAAegWAQ4AMUiBzcuagecDjHZfkXNnzfLTyWT7ecPVXNwhw4lA4JrcnDrs4G+TXNwycR6DzvZPvldMqm+bFzBAShWJwHO9rtsP2D7fts3217dRT0AlK31AGf7WZLeLmk2Il4saVrSZW3XA0D5usrBrZC0xvYxSWsl/aSjeoy/pmO8kryca8aNpZPwk0nw2Ti4LAc3dzzPwR05Un3NmKuWuenE+kEm26c7j5FvW25av4KLiB9L+qCkH0raL+l/I+L2tusBoHxd3KJuknSppHMkPVPSOttvTo7baXu37d2PHnis7WoCKEAXnQy/J+l7EfFoRByTdJuk31x8UETsiojZiJjdumVz65UEMPm6CHA/lPQy22vdG4R0kaS9HdQDQOFa72SIiLtt3yrpXklzkr4uaVfb9ZgYjTsZkgT89Ez+mtmA2ah2FGSdBE8kZceT50rSkaNJJ0O2A9dUtoNWUve6CfSNB/piuemkFzUi3ifpfV2cG8DycdJbVNsEIgATqckV3Ptsr5V0mnq3lbdExP+MtloAMLwmAS4kHZb0BUkXSvoP25dFxJ6R1myZqZv0nS6Eme52n+TVZlbmJ6ubhL/I0WPVfNmhJIc2pwEG+mYT85O6O6v79ACLB2Q7j5GXW3aavNMf7OfMJOlW2zdI+qikV4+sVgCwBJoMEzlg+yVPfhMR35G0dXRVAoCl0eQK7u2SbrF9j6RvSjpP0vdGWisAWAInDXARscf2+erNQHixpC9JunnUFUPfEJPt63NwNePjFplLcnBPJAtjztWMgzt2NBnzlh2b1Ser+7CT7bHsNMo2R8QRSf/c/wKAicCKvgCKRYADUCwCHIBisavWJGo60HdFTSdDlsBPXvPYXLWT4OAAnQzZQOG07jNJJ0NW97QzQWL1XtThCg5AsQhwAIpFgANQLHJwYy6bIB7p5PJkV62ayemxclW1MJmAfzzZLevxZLL9wWwRS0nHknJni1uuSrbFXZmUDTDQl4n1kLiCA1AwAhyAYhHgABSLHNwkajoOrm7TmVVrqmXJWLT5JAf382Rs3KFkbJyUb0bjFUk9k/o4m4Bfl4NjYj1qdHIFZ/tU27faftD2Xtsv76IeAMrW1RXcRyT9S0T8oe2VktZ2VA8ABWs9wNneKOlVkv5YkiLiqKSjbdcDQPm6uEV9tqRHJf297a/b/pjtdR3UA0DhurhFXaHe7lxv6+9y/xFJV0n6y4UH2d4paackbd+2rfVKTpxs8G/dyr2rk8+T1cnA2sQTyeDdx5POCEmq7qkleVUyiX5NkqHIOkjqOhkY1IsaXVzB7ZO0LyLu7n9/q3oB7ykiYldEzEbE7NYtm1utIIAytB7gIuK/Jf3I9rn9ooskfavtegAoX1e9qG+TdFO/B/W7kv6ko3oAKFgnAS4i7pM028W5S5BOwM8Wg6zJWXnt+urz11QH205NV89zMNmZPlsEU6q5PUjOozUbqmVZ/jDLM4qJ9ajHVC0AxSLAASgWAQ5AsQhwAIrFaiLFSBLtNSv6Zkl9r6+WrVpZ7aTIdtA6VrOr1nSW/F9XHWTsNdVOj7TuNZ0MQB3eMQCKRYADUCwCHIBikYMrRZbvqsnBOVlBNzaeWilbs6aag5saYPXc1TPJQONTTqmWrdtYLUtzcAzoxWC4ggNQLAIcgGIR4AAUixxcIfIJ+DULRK5IdrbfVF1zb93G6iKY65MJ+Ifn88/JjRuqE+adnMerkgUvp6pvTSbVY1BcwQEoFgEOQLEIcACKRYADUCw6GYpWk5SfSXa2Ou2MStHKZ1YH/26d2Vcpm6o5zaZNSWfGlmdUy1YmO3plKxQDA+JLine truncated
|
||||||
|
"text/plain": [
|
||||||
|
"<Figure size 432x288 with 1 Axes>"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"metadata": {
|
||||||
|
"needs_background": "light"
|
||||||
|
},
|
||||||
|
"output_type": "display_data"
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"likelihood = stats.norm(mm, ss).pdf(lengths[0]) * stats.norm(mm, ss).pdf(lengths[1])\n",
|
||||||
|
"aspect = mm.max() / ss.max() / 3\n",
|
||||||
|
"extent = [mm.min(), mm.max(), ss.min(), ss.max()]\n",
|
||||||
|
"# extent = left right bottom top\n",
|
||||||
|
"plt.imshow(likelihood, cmap='Reds', aspect=aspect, extent=extent) \n",
|
||||||
|
"plt.xlabel(r'$\\mu$')\n",
|
||||||
|
"plt.ylabel(r'$\\sigma$')"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 14,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"text/plain": [
|
||||||
|
"Text(0, 0.5, '$\\\\sigma$')"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"execution_count": 14,
|
||||||
|
"metadata": {},
|
||||||
|
"output_type": "execute_result"
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAATgAAAEMCAYAAABQlKrRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAG11JREFUeJzt3W+sZVV5BvDnuffOgCBUcAalMjDUKok1FOHW2GptK9WiNeKHfoBoQyvJpCb1X2oVSlqTflIxVZM2MROlYiSQFjGaVivUSomt0g7ICDggVPyDUpmJTdU0zNw/bz+cM3K551lz3z377r3PWff5JTdz77r7nL3WPmfW3edd71qLEQEzsxrNDV0BM7OuuIMzs2q5gzOzarmDM7NquYMzs2q5gzOzarmDM7NquYMzs2q5gzOzai0MXYGMHTt2xDm7z97U5+xmAkfuSfVR7UrVjJQoPTp5bMRq6rjVwsXUz6mOVfXJY7KUFGWFR88ljyUn7xH0caVaZo9V9SlpcuzxH5l+xs1/Stx919cORcTOjY6biQ7unN1n49/u/PKmPmeTDkE+Pt3JiE5CHLcaK/I8K+Lx6tiV1eWJsuWYLAOAIyuHJ8tWj4gycdyKOm6yrHTskqjniqhnkymEqkOY5+Rbe9vcZNn2+e3yObfPTZarY7fPnZB87ORxALAg6jkv6jnH+cnjCh/A1LG6cxedc4M/AkqTzr2Npy2c/J3Mcf6IambVcgdnZtWaiY+otWnyUVg/QfbjdT42tqribeqjtPzIPFlWKpcfr5PnLpEfizB5ntWY/HterHs6tJCLUzZ5LWSAuMGnPPWcTT5m1sR3cGZWLXdwZlYtd3BmVi13cGZWLQ8yTIliSF0OKCSfs0kCbouBh5VSDp8cZMgNPJSC/8qc+Du9KoLqqj6lurcZUMgnOEMOHsgj1fugMG6wNYcTNN/BmVm1OuvgSF5H8nGS94nfvZNkkNzR1fnNzLq8g/s4gEvWF5LcBeCVAL7b4bnNzLqLwUXEHSR3i199EMC7AHymq3NPky7m9GcThVWyKtBgsn0yLldO9M0l9S6vTh7XJBlaTfbn3OTf7jkZ65uctzkq39x4W+m1UFongsvnnLQVYnW9xuBIvg7A9yNif5/nNbOtqbdRVJInAbgGwKuSx+8BsAcAdp29q8OamVmt+ryDey6AcwHsJ/ltAGcBuJvks9XBEbE3IhYjYnHnTo9FmFlzvd3BRcS9AM44+vO4k1uMiEN91WGaFRenTC+imc+90hPJszlrKpcsP9lerf2WzY0rUeufqfPMi8n2Tequ8wLFccmJ+kDL17eTCG9dukwTuRHAVwCcR/JRkld2dS4zM6XLUdTLN/j97q7ObWYGeCaDmVXMHZyZVcuT7Y9Xcneojk4+WdLBrloqWK4n0OsBAbXpjUz0bbnpzKpawVbMRFfnmS8k+sp2qsGDDl4LcLj3kc4Int2UYN/BmVm13MGZWbXcwZlZtRyDG0CTXZf043PP2SS5NJvY2mRxSpnoKzZ+bjKBX1ELXqrzzItd6JssFJBO/m35WrR7J5R21dqafAdnZtVyB2dm1XIHZ2bVcgdnZtXaEoMMTZJGOzh5u4fLZMzcIEUpuTQ/oNBuNRG1Uu+yCN4vJwceSubE4IH60z0n6rMw1241kTYDD0d/M1mkdtBq+R6WCcX9DD2o/3/s6dy+gzOzarmDM7NquYMzs2ptiRjcLGi2s33LCd4tdtCSE+hFDG30+FyyrNpdfmU1H4Obn5us+5xYvVfu8lWo+/JcLi6Y3mmrycIHKt7WYGd7e5Lv4MysWu7gzKxa7uDMrFqOwW2gk52LOsjLy8bQRqfP5rdlY2iFPLhkztuSyoNrEINbFX+nKaaXq3y5hUY7guWuR/b6Aj2+v1rE60p1VNd42nS5q9Z1JB8ned+asmtJPkDy6yQ/TfIZXZ3fzKzLj6gfB3DJurLbALwwIs4H8E0AV3d4fjPb4jrr4CLiDgA/Wld2a8TP8gy+itHu9mZmnRhykOFNAD4/4PnNrHKDDDKQvAbAMoAbjnHMHgB7AGDX2bt6qlk7TRI8Jw8sJeVu/rnbrOirBhRU8B3QAwpyAr4oa5ToK9qjVvmdo0hSntN117uHtVjRt8FroY9ThYXHith/tk6zMHDQRO93cCSvAPBaAG+IYyzzERF7I2IxIhZ37tzRXwXNrBq93sGRvATAuwH8RkT8X5/nNrOtp8s0kRsBfAXAeSQfJXklgL8GcAqA20jeQ/IjXZ3fzKyzO7iIuFwUf6yr8/WtTXpmk8e2mRjfZLerdLxNxtX0hHUVm9NxOTWBPx+DCxVvE7Gk+bnJ44p1V3Xa5EUwS+WtYrkobE7fw2OnkadqmVm13MGZWbXcwZlZtTzZfko02ZRELpKYLAPa7WKfjUMBehL9Ugd5cDqfS8TlRH22Fzad0e1ssdt9Kc8x+7qJRTA7mahfGd/BmVm13MGZWbXcwZlZtdzBmVm1tuwgQxcB2jYT45s8Zzr5F01Wq1VJvWqQQe1WX9iZSjx+aXVpouzIymRZaZVgZV7soCWPWxWJvoW6652+jn/QpfRa5Hflymszgb9tUu+0TeD3HZyZVcsdnJlVyx2cmVVry8bgZlsy6bO4iKaarJ+LL2UXsSwdq5J/9SKYOjamxNzk25icjPssrM6nzl06fzb5t+1r0W4pB1vLd3BmVi13cGZWLXdwZlYtx+A6lo2xFHOXkpOxm+Re5Re8FDvOJ+Nyo3KRB7cyeWzbPLjV8tYeT7FNxOpKdW9zPdq+FnKyvXofiAn449+IEpWfVj/fwZlZtdzBmVm13MGZWbW63FXrOpKPk7xvTdnpJG8j+dD439O6Or+ZWZeDDB/HaJvAT6wpuwrAFyPivSSvGv/87g7r0F4xgL3JyZhNknJFwFpO0G6woq+adC4TfcVxR8TAwah8cvBABfWPiLLVQgKusjI3mcA7JxJ9VX22zW2Tz7m8cPzXo9GKvi1e3yY72+ep8xQOFdd42nR2BxcRdwD40briSwFcP/7+egCv7+r8ZmZ9x+CeFRGPAcD43zN6Pr+ZbSFTO8hAcg/JfST3HTx4aOjqmNkM6jvR94ckz4yIx0ieCeDx0oERsRfAXgC4aPHC2Zh9nEw4bduY7E5MrXe2bzkxXh2r4nVLMtE3f5XmRd3VzvYLDRJ99U5fLRJ9C69Fk93QstKPlrt3TX9crYm+7+A+C+CK8fdXAPhMz+c3sy2kyzSRGwF8BcB5JB8leSWA9wJ4JcmHALxy/LOZWSc6+4gaEZcXfnVxV+c0M1vLk+27lp04XXp4doOZBhO89UKWajOZXLyttGikyjs7IjadObzSbmf7+bnJDyLk5HkWVibz5U6Y3y6fU7ZTXo/kRjRNJtu33YhGxvBUvK30BPWY2lFUM7O23MGZWbXcwZlZtdzBmVm1PMgwJUor+mYTQbMDD0BpZdrJYLl6vErUVQMHo3I1MV8l/7YcZBDXQ022X5hrUHfRTp3AqwYZco8F8oskyLLiir52lO/gzKxa7uDMrFru4MysWlXF4NpOUm50rhbJmE12M0+fp8EEb71wYzapt93O9tkY3JIoK9k2n1vwcvt8u7pnk39lom9xsv3xJ/o2e89kj1O7b21+RnDp/yo3ebK/7+DMrFru4MysWu7gzKxa7uDMrFpVDTK0VUq27efc+d9kVw4ptSeb1KtX71WDBKVk2cljDy9PBurVaiJNBhlU+F6FqreLwYhS3bNtzyb/Fl8LUXs9INFgNZGe9DUg0Ybv4MysWu7gzKxa7uDMrFqOwW2g9Q5YLZ8hO7G+7a5a2cRWNTld7YoF5JN6nxBxuSaT7dUOXOovt04oLi0UMFmevUatd9VKvr5NtHkfFje2P+5n7I/v4MysWoN0cCTfQfJ+kveRvJHkiUPUw8zq1nsHR/I5AN4KYDEiXghgHsBlfdfDzOo3VAxuAcDTONr66CQAPxioHoOQ8ZBCjKXVrlqFuI/eCUrEy2QcSh2nc9ZUbE2V6d3u8zG4bfOTf6fVX241Kf/EBf1fQLdz8nqo65bNlyuVpyfbl94zXgjzZ3q/g4uI7wP4AIDvAngMwP9GxK1918PM6jfER9TTAFwK4FwAPw/gZJJvFMftIbmP5L6DBw/1XU0zq8AQgwy/DeCRiDgYEUsAbgHwa+sPioi9EbEYEYs7d+7ovZJmNvuG6OC+C+AlJE/iaHW7iwEcGKAeZla53gcZIuJOkjcDuBvAMoCvAdjbdz3ayodxc0cWd9VqMdleBcBHj88NMmQn1qsJ9ICeRK/KnlgW524w2X4lJgcP1KTv7cnJ/4BuZ3ZwJnt9gfxk+/wqvyXZ9+GkWUjoLRlkFDUi3gPgPUOc28y2jg0/opJ0R2RmMylzB/cekicBOB2jj5U3RcT/dFstM7P2Mh1cAHgCwBcAXAjg30leFhH7O63ZJuprIcv0eRpMnE7vci7KVko7RiXjbUvZCfiFeNmRbKLvkphs3yDRd3VVLbw46QS14GUx0VckHyd3FJM7bRVeizavb5E6VlwQvWDl5htyYcxMB/fAOGYGADeT/DiAjwB4RWe1MjPbBJk0kUMkLzr6Q0R8E8DO7qpkZrY5MndwbwVwE8m7ANwL4HwAj3RaKzOzTbBhBxcR+0legNEMhBcC+BKAG7uuWC2y+5GXN4jJ5bzpPLi2C16qPLjcpjEA8IScRC9iViIPbkWUlYTIg5ubm4zxqPqU6q7aqTeiyS14WXwtkq9lNjduVJ4r2wpSeXARcRjAP46/zMxmglf0NbNquYMzs2q5gzOzanlXrQ0VwrMtdzmafLr8ZHt17uwO60B+QEHtOJXdKQsADouBgqUlEagXib5q4KFEXbvRQjXr6jOfr3t2B678wENp4QMx+NB2p602ObQySXh299XyHZyZVcsdnJlVyx2cmVXLMbi12sbV1OOz8ZQmC14mE0Gb7OSkJpIvifiSSowtxbFkUq+Ity31FINLine truncated
|
||||||
|
"text/plain": [
|
||||||
|
"<Figure size 432x288 with 1 Axes>"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"metadata": {
|
||||||
|
"needs_background": "light"
|
||||||
|
},
|
||||||
|
"output_type": "display_data"
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"prior = stats.norm(200, 15).pdf(mm) * stats.cauchy(0, 10).pdf(ss)\n",
|
||||||
|
"\n",
|
||||||
|
"plt.imshow(prior, cmap='Greens', aspect=aspect, extent=extent)\n",
|
||||||
|
"plt.xlabel(r'$\\mu$')\n",
|
||||||
|
"plt.ylabel(r'$\\sigma$')"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 15,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"text/plain": [
|
||||||
|
"Text(0, 0.5, '$\\\\sigma$')"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"execution_count": 15,
|
||||||
|
"metadata": {},
|
||||||
|
"output_type": "execute_result"
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAATgAAAEMCAYAAABQlKrRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFGxJREFUeJzt3WuMXVd5xvHnmfHdsXMbh5jEjgOFSDRKAx4QBUovKTRQRPjQD4mgSkskq0jlplKaNGqR+gkIKiC1ErIgTSqiRG0IArVQklLSqC24dUxMHBxIytXgYltQ7MS3ubz9cHbKZPa77X2wz9nnrPn/pCOfWbPPWWvZ43fWeddlOyIEACWa6LoBADAoBDgAxSLAASgWAQ5AsQhwAIpFgANQLAIcgGIR4AAUiwAHoFjLum5AG1NTU3HZZVu6bsZADWI/SdMmlWhZW9tNLme8Gcatippfnr6+3Zv2U09bg3hPPNuuXQ8fiogNp7tuLALcZZdt0b/v2Nl1MwaqactcVjyfFGbXzTW859x8vbxt2XxW1hRIs/qTaDSRRISJ5LrJ7EJJy5LyiaQse33DW8ot29n2tTi7Vi/3d9tcx0dUAMUiwAEo1lh8RMXPqY/cWPZpMv94nF3X8PG6ZUXzSdaqn095nIeDJozgABSLAAegWAQ4AMUiwAEoFpMMhegr0d52AW/LssZ1cC2bM5lc2XbSo6fdjETbdXkoByM4AMUaWICzfbvtA7b3JN97j+2wPTWo+gFgkCO4OyRdu7jQ9iZJr5H0vQHWDQCDy8FFxEO2tyTf+rCk90r6zKDqRk9TyirNraX7W9uVZXtje9fWy7KUV7rQN8vBpbU0LTQeTm6Nfaejbag5ONtvlPSDiNg9zHoBLE1Dm0W1vUbSrZJe2/L6bZK2SdKmzZsH2DIApRrmCO75ki6XtNv2dyRdKmmX7YuziyNie0RMR8T0hqnTHvsEADVDG8FFxKOSLnrm6yrITUfEoWG1AT1pbi25Llvf1rasqZ70IMqk9onkun4O1jzTNW9k1sowyGUid0v6sqQrbO+zfdOg6gKAzCBnUW84zfe3DKpuAJDYyQCgYAQ4AMVis33B+jlpN50oyBb1Jhf2c8OcrPZIJxTaXZe/I9DDCA5AsQhwAIpFgANQLHJwhWjKg+XXtivMF/W2u0F0c+VJWfJrNiLJtzX2sd2yXDbGLz2M4AAUiwAHoFgEOADFIsABKBaTDAXrZwFs2wmFMz1NZCJJ9Kd30MrfsqGePi7GksIIDkCxCHAAikWAA1AscnBL0Jmc6JveVashCde2NL2DVh+b7YEmjOAAFIsAB6BYBDgAxSIHV4iWZ0s2Fre9i322sb6fO9tne+iztXF5/i+tBmg0yLtq3W77gO09C8pus/247a/Z/rTt8wZVPwAM8iPqHZKuXVT2gKQrI+IqSd+UdMsA6wewxA0swEXEQ5J+vKjs/oiYrb78inp3tweAgehykuGtkj7fYf0ACtfJJIPtWyXNSrrrFNdsk7RNkjZt3jyklpWlKSefJeuzpH42odC2rKmeiYlsQqHdBEfXcwycCDx+hj6Cs32jpDdIenOc4pztiNgeEdMRMb1hasPwGgigGEMdwdm+VtKfSPrViDg6zLoBLD2DXCZyt6QvS7rC9j7bN0n6K0nrJD1g+xHbHxtU/QAwsBFcRNyQFH9iUPWhvbYLfdsebtlPDi4zn+Tl2t7560w15dVIt5WBrVoAikWAA1AsAhyAYrHZvmBNq3Ba59sGsA4uM5m8/kyzbVlujbTa0sMIDkCxCHAAikWAA1AsAhyAYjHJsASlJ+2mEw/1stm5+VrZzFw+JTCfTBWsSH6njstme4wfRnAAikWAA1AsAhyAYpGDK0WWV2t/6Rkdbnlifi6vJzvwMlluu7zlHbTOdK99toGexb9lYwQHoFgEOADFIsABKBY5uII15azyu9jXr8vybSfn6+vgTszWy5pMJomwlfPt1sb1g9waJEZwAApGgANQLAIcgGIN8q5at9s+YHvPgrILbD9g+4nqz/MHVT8ADHIEd4ekaxeV3SzpixHxAklfrL6GeifQZo9BiKg/5iNqj7n5+uP47Hz9MTfX+nFibr72yOru2rD+LTBYAwtwEfGQpB8vKr5O0p3V8zslvWlQ9QPAsHNwz4mI/ZJU/XnRkOsHsISM7CSD7W22d9reefDQwa6bA2AMDXuh749sb4yI/bY3SjrQdGFEbJe0XZK2bp3uPilTkLab7bODLE/M1TfWZ2VNliV3sZ+br/8YnnEarmXKjNRa2YY9gvuspBur5zdK+syQ6wewhAxymcjdkr4s6Qrb+2zfJOn9kl5j+wlJr6m+BoCBGNhH1Ii4oeFb1wyqTgBYiM32Iy7NEWX5qT5ySW3vbD+TbKw/nuTbnp7Jc3BJuk3LJ+ofGmZb3vCmCQdZosnIzqICwJkiwAEoFgEOQLEIcACKxSTDEpTl7+ezu2Uld7E/PptMMpzMT/TNJhlWLatfmy0y7kfbCQU2zC89jOAAFIsAB6BYBDgAxSIHtwRlGa9sYW22if7IyXrZ4eP5Qt/lk/Wc1zkr6zm42Wh3V66mFFqWW0vLWtWCkjCCA1AsAhyAYhHgABSLAAegWEwyjKF8wWr7xbLZaSLZYtuTyWkiT52ol/30+Gxaz8pl9d+f61dNJnWnL69pmiTIyln8C4kRHICCEeAAFIsAB6BY5OAKcaaZpGy/e7ax/siJpKxhoe/xZfU3vXBNPeHW+vTehnxZeqJvyzKUjREcgGJ1EuBsv9v2Y7b32L7b9qou2gGgbEMPcLYvkfQOSdMRcaWkSUnXD7sdAMrXVQ5umaTVtmckrZH0w47aMfLSTeNJzupM13Nl6+COJjm4/z1WX/N2+NjJ9D1XLa+veTs6s7xed8scXHaApsRaNjQb+gguIn4g6UOSvidpv6SfRsT9w24HgPJ18RH1fEnXSbpc0nMlrbX9luS6bbZ32t558NDBYTcTQAG6mGT4TUnfjoiDETEj6T5Jr1h8UURsj4jpiJjeMLVh6I0EMP66CHDfk/Ry22vcS55cI2lvB+0AULihTzJExA7b90raJWlW0lclbR92O5ayLKefJfqPJXfLeirZWH/46Exaz+zK+nsen62/ZyQHBbQ9pbdXnpUx8YCOZlEj4n2S3tdF3QCWjtN+RLVNIAIwltqM4N5ne42kC9T7WHlPRPxksM0CgDPXJsCFpOOSviDpJZL+w/b1EbF7oC1DX9JDH/tIQ2V3tjo6Uy87fLS+qPfw0/lC32wT/fHkPTPZot7Ghb5ty8jLLTltAtzjVc5Mku61fYekj0n6jYG1CgDOgjbLRA7Z3vrMFxHxTUksTAMw8tqM4N4h6R7bD0t6VNJVkr490FYBwFlw2gAXEbttX63eDoQrJX1J0t2DbhiatT3McaKPdWNZvuzpZB3c4WP1NW9HjpxoXc/JubYb6+svbupPUznQah1cRJyQ9I/VAwDGAif6AigWAQ5AsQhwAIrFXbUK1k/qPduAny30fTpZ1Ns0ybAsubP9bDLJMJlNKCSreicbVvpyBy00YQQHoFgEOADFIsABKBY5uDGUHwaZXZe/vu3C2GMn63fVeuqpeg7u6SNH09evWlX/8cru3pW1J8u3NTW7n8MxsbQwggNQLAIcgGIR4AAUixxcIbKMUz+b7TMzc+3WwR1tyMGtXbemVT0rJuu/Z5dl6+Ca+tOqFixFnYzgbJ9n+17bj9vea/uXu2gHgLJ1NYL7qKR/iojfsb1CUrtf9QDQh6EHONvrJb1a0u9JUkSclJQf6g8AZ6CLj6jPk3RQ0t/Y/qrtj9te20E7ABSui4+oy9S7O9fbq7vcf1TSzZL+bOFFtrdJ2iZJmzZvHnoji9CQfc8W0WYTEjPJXeiPPX2sft2RI2k9J0+cW29S0qaVE5O1smWT7TbgN70nIHUzgtsnaV9E7Ki+vle9gPcsEbE9IqYjYnrDFPe4AdC/oQe4iPgfSd+3fUVVdI2krw+7HQDK19Us6tsl3VXNoH5L0u931A4ABeskwEXEI5Kmu6i7VP3cVSvLwWULa+eSUzBPHEsOt3z6J2k9c7Mba2UrkkMws7KsPU13tgeasFULQLEIcACKRYADUCwCHIBicZpIIbITbCec3CpL+STDiol2v+tmTszUC48dTq+dn68vFF67ol7PymySITlhpOmUXk7vRRNGcACKRYADUCwCHIBikYMrWFNuKsvBLU9yXtlJu3MzSQ7uZH0DviRNTtY30Z+3uv4jt3I5C30xGIzgABSLAAegWAQ4AMUiB1ewppRVlt/K8mCrV7T88Ug25UvS8pXLa2Ub1mY5uHquLssTAv1iBAegWAQ4AMUiwAEoFgEOQLGYZChY0x70dLN9suH9/DX1H49lK1bUymaX1cskaf36VbWyqdUra2XZZvvsDlpsqke/GMEBKBYBDkCxOgtwtierO9v/Q1dtAFC2LnNw75S0V9L6DttQtH42269KFttevK6+UHfNujW1ssPrptJ6NmxYWyu7aHU9L5fm4Ei34SzoZARn+1JJvy3p413UD2Bp6Ooj6kckvVdS/UxrADhLhh7gbL9B0oGIePg0122zvdP2zoOHDg6pdQBK0kUO7pWS3mj79ZJWSVpv+5MR8ZaFF0XEdknbJWnr1ul8Nzd+Ltkd71clm+0vWVfPl12woZ4yPX50U1rP8y9eVytbv6ae11ue5OBY84azYegjuIi4JSIujYgtkq6X9C+LgxsAnA2sgwNQrE63akXEg5Ie7LINAMrFCA5AsdhsvwRl+fvsVN1LzlldK9uy+bxa2WRy9y1JuvqSc2pl61bVf+Q4vReDwggOQLEIcACKRYADUCxycEtQtoh2+WS97Ly19YMsX/q882tlU8mCYEm6cqq+KHjNynqujxQcBoURHIBiEeAAFIsAB6BY5OAgKV+Ldu7q+o/HKy6t5+BeOFU/BFOSNp7b7nBLNtZjUBjBASgLine truncated
|
||||||
|
"text/plain": [
|
||||||
|
"<Figure size 432x288 with 1 Axes>"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"metadata": {
|
||||||
|
"needs_background": "light"
|
||||||
|
},
|
||||||
|
"output_type": "display_data"
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"unnormalized_posterior = prior * likelihood\n",
|
||||||
|
"posterior = unnormalized_posterior / np.nan_to_num(unnormalized_posterior).sum()\n",
|
||||||
|
"plt.imshow(posterior, cmap='Blues', aspect=aspect, extent=extent)\n",
|
||||||
|
"plt.xlabel(r'$\\mu$')\n",
|
||||||
|
"plt.ylabel(r'$\\sigma$')"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 7,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"text/plain": [
|
||||||
|
"Text(0, 0.5, '$\\\\sigma$')"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"execution_count": 7,
|
||||||
|
"metadata": {},
|
||||||
|
"output_type": "execute_result"
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAATgAAAEMCAYAAABQlKrRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFGxJREFUeJzt3WuMXVd5xvHnmfHdsXMbh5jEjgOFSDRKAx4QBUovKTRQRPjQD4mgSkskq0jlplKaNGqR+gkIKiC1ErIgTSqiRG0IArVQklLSqC24dUxMHBxIytXgYltQ7MS3ubz9cHbKZPa77X2wz9nnrPn/pCOfWbPPWWvZ43fWeddlOyIEACWa6LoBADAoBDgAxSLAASgWAQ5AsQhwAIpFgANQLAIcgGIR4AAUiwAHoFjLum5AG1NTU3HZZVu6bsZADWI/SdMmlWhZW9tNLme8Gcatippfnr6+3Zv2U09bg3hPPNuuXQ8fiogNp7tuLALcZZdt0b/v2Nl1MwaqactcVjyfFGbXzTW859x8vbxt2XxW1hRIs/qTaDSRRISJ5LrJ7EJJy5LyiaQse33DW8ot29n2tTi7Vi/3d9tcx0dUAMUiwAEo1lh8RMXPqY/cWPZpMv94nF3X8PG6ZUXzSdaqn095nIeDJozgABSLAAegWAQ4AMUiwAEoFpMMhegr0d52AW/LssZ1cC2bM5lc2XbSo6fdjETbdXkoByM4AMUaWICzfbvtA7b3JN97j+2wPTWo+gFgkCO4OyRdu7jQ9iZJr5H0vQHWDQCDy8FFxEO2tyTf+rCk90r6zKDqRk9TyirNraX7W9uVZXtje9fWy7KUV7rQN8vBpbU0LTQeTm6Nfaejbag5ONtvlPSDiNg9zHoBLE1Dm0W1vUbSrZJe2/L6bZK2SdKmzZsH2DIApRrmCO75ki6XtNv2dyRdKmmX7YuziyNie0RMR8T0hqnTHvsEADVDG8FFxKOSLnrm6yrITUfEoWG1AT1pbi25Llvf1rasqZ70IMqk9onkun4O1jzTNW9k1sowyGUid0v6sqQrbO+zfdOg6gKAzCBnUW84zfe3DKpuAJDYyQCgYAQ4AMVis33B+jlpN50oyBb1Jhf2c8OcrPZIJxTaXZe/I9DDCA5AsQhwAIpFgANQLHJwhWjKg+XXtivMF/W2u0F0c+VJWfJrNiLJtzX2sd2yXDbGLz2M4AAUiwAHoFgEOADFIsABKBaTDAXrZwFs2wmFMz1NZCJJ9Kd30MrfsqGePi7GksIIDkCxCHAAikWAA1AscnBL0Jmc6JveVashCde2NL2DVh+b7YEmjOAAFIsAB6BYBDgAxSIHV4iWZ0s2Fre9i322sb6fO9tne+iztXF5/i+tBmg0yLtq3W77gO09C8pus/247a/Z/rTt8wZVPwAM8iPqHZKuXVT2gKQrI+IqSd+UdMsA6wewxA0swEXEQ5J+vKjs/oiYrb78inp3tweAgehykuGtkj7fYf0ACtfJJIPtWyXNSrrrFNdsk7RNkjZt3jyklpWlKSefJeuzpH42odC2rKmeiYlsQqHdBEfXcwycCDx+hj6Cs32jpDdIenOc4pztiNgeEdMRMb1hasPwGgigGEMdwdm+VtKfSPrViDg6zLoBLD2DXCZyt6QvS7rC9j7bN0n6K0nrJD1g+xHbHxtU/QAwsBFcRNyQFH9iUPWhvbYLfdsebtlPDi4zn+Tl2t7560w15dVIt5WBrVoAikWAA1AsAhyAYrHZvmBNq3Ba59sGsA4uM5m8/kyzbVlujbTa0sMIDkCxCHAAikWAA1AsAhyAYjHJsASlJ+2mEw/1stm5+VrZzFw+JTCfTBWsSH6njstme4wfRnAAikWAA1AsAhyAYpGDK0WWV2t/6Rkdbnlifi6vJzvwMlluu7zlHbTOdK99toGexb9lYwQHoFgEOADFIsABKBY5uII15azyu9jXr8vybSfn6+vgTszWy5pMJomwlfPt1sb1g9waJEZwAApGgANQLAIcgGIN8q5at9s+YHvPgrILbD9g+4nqz/MHVT8ADHIEd4ekaxeV3SzpixHxAklfrL6GeifQZo9BiKg/5iNqj7n5+uP47Hz9MTfX+nFibr72yOru2rD+LTBYAwtwEfGQpB8vKr5O0p3V8zslvWlQ9QPAsHNwz4mI/ZJU/XnRkOsHsISM7CSD7W22d9reefDQwa6bA2AMDXuh749sb4yI/bY3SjrQdGFEbJe0XZK2bp3uPilTkLab7bODLE/M1TfWZ2VNliV3sZ+br/8YnnEarmXKjNRa2YY9gvuspBur5zdK+syQ6wewhAxymcjdkr4s6Qrb+2zfJOn9kl5j+wlJr6m+BoCBGNhH1Ii4oeFb1wyqTgBYiM32Iy7NEWX5qT5ySW3vbD+TbKw/nuTbnp7Jc3BJuk3LJ+ofGmZb3vCmCQdZosnIzqICwJkiwAEoFgEOQLEIcACKxSTDEpTl7+ezu2Uld7E/PptMMpzMT/TNJhlWLatfmy0y7kfbCQU2zC89jOAAFIsAB6BYBDgAxSIHtwRlGa9sYW22if7IyXrZ4eP5Qt/lk/Wc1zkr6zm42Wh3V66mFFqWW0vLWtWCkjCCA1AsAhyAYhHgABSLAAegWEwyjKF8wWr7xbLZaSLZYtuTyWkiT52ol/30+Gxaz8pl9d+f61dNJnWnL69pmiTIyln8C4kRHICCEeAAFIsAB6BY5OAKcaaZpGy/e7ax/siJpKxhoe/xZfU3vXBNPeHW+vTehnxZeqJvyzKUjREcgGJ1EuBsv9v2Y7b32L7b9qou2gGgbEMPcLYvkfQOSdMRcaWkSUnXD7sdAMrXVQ5umaTVtmckrZH0w47aMfLSTeNJzupM13Nl6+COJjm4/z1WX/N2+NjJ9D1XLa+veTs6s7xed8scXHaApsRaNjQb+gguIn4g6UOSvidpv6SfRsT9w24HgPJ18RH1fEnXSbpc0nMlrbX9luS6bbZ32t558NDBYTcTQAG6mGT4TUnfjoiDETEj6T5Jr1h8UURsj4jpiJjeMLVh6I0EMP66CHDfk/Ry22vcS55cI2lvB+0AULihTzJExA7b90raJWlW0lclbR92O5ayLKefJfqPJXfLeirZWH/46Exaz+zK+nsen62/ZyQHBbQ9pbdXnpUx8YCOZlEj4n2S3tdF3QCWjtN+RLVNIAIwltqM4N5ne42kC9T7WHlPRPxksM0CgDPXJsCFpOOSviDpJZL+w/b1EbF7oC1DX9JDH/tIQ2V3tjo6Uy87fLS+qPfw0/lC32wT/fHkPTPZot7Ghb5ty8jLLTltAtzjVc5Mku61fYekj0n6jYG1CgDOgjbLRA7Z3vrMFxHxTUksTAMw8tqM4N4h6R7bD0t6VNJVkr490FYBwFlw2gAXEbttX63eDoQrJX1J0t2DbhiatT3McaKPdWNZvuzpZB3c4WP1NW9HjpxoXc/JubYb6+svbupPUznQah1cRJyQ9I/VAwDGAif6AigWAQ5AsQhwAIrFXbUK1k/qPduAny30fTpZ1Ns0ybAsubP9bDLJMJlNKCSreicbVvpyBy00YQQHoFgEOADFIsABKBY5uDGUHwaZXZe/vu3C2GMn63fVeuqpeg7u6SNH09evWlX/8cru3pW1J8u3NTW7n8MxsbQwggNQLAIcgGIR4AAUixxcIbKMUz+b7TMzc+3WwR1tyMGtXbemVT0rJuu/Z5dl6+Ca+tOqFixFnYzgbJ9n+17bj9vea/uXu2gHgLJ1NYL7qKR/iojfsb1CUrtf9QDQh6EHONvrJb1a0u9JUkSclJQf6g8AZ6CLj6jPk3RQ0t/Y/qrtj9te20E7ABSui4+oy9S7O9fbq7vcf1TSzZL+bOFFtrdJ2iZJmzZvHnoji9CQfc8W0WYTEjPJXeiPPX2sft2RI2k9J0+cW29S0qaVE5O1smWT7TbgN70nIHUzgtsnaV9E7Ki+vle9gPcsEbE9IqYjYnrDFPe4AdC/oQe4iPgfSd+3fUVVdI2krw+7HQDK19Us6tsl3VXNoH5L0u931A4ABeskwEXEI5Kmu6i7VP3cVSvLwWULa+eSUzBPHEsOt3z6J2k9c7Mba2UrkkMws7KsPU13tgeasFULQLEIcACKRYADUCwCHIBicZpIIbITbCec3CpL+STDiol2v+tmTszUC48dTq+dn68vFF67ol7PymySITlhpOmUXk7vRRNGcACKRYADUCwCHIBikYMrWFNuKsvBLU9yXtlJu3MzSQ7uZH0DviRNTtY30Z+3uv4jt3I5C30xGIzgABSLAAegWAQ4AMUiB1ewppRVlt/K8mCrV7T88Ug25UvS8pXLa2Ub1mY5uHquLssTAv1iBAegWAQ4AMUiwAEoFgEOQLGYZChY0x70dLN9suH9/DX1H49lK1bUymaX1cskaf36VbWyqdUra2XZZvvsDlpsqke/GMEBKBYBDkCxOgtwtierO9v/Q1dtAFC2LnNw75S0V9L6DttQtH42269KFttevK6+UHfNujW1ssPrptJ6NmxYWyu7aHU9L5fm4Ei34SzoZARn+1JJvy3p413UD2Bp6Ooj6kckvVdS/UxrADhLhh7gbL9B0oGIePg0122zvdP2zoOHDg6pdQBK0kUO7pWS3mj79ZJWSVpv+5MR8ZaFF0XEdknbJWnr1ul8Nzd+Ltkd71clm+0vWVfPl12woZ4yPX50U1rP8y9eVytbv6ae11ue5OBY84azYegjuIi4JSIujYgtkq6X9C+LgxsAnA2sgwNQrE63akXEg5Ie7LINAMrFCA5AsdhsvwRl+fvsVN1LzlldK9uy+bxa2WRy9y1JuvqSc2pl61bVf+Q4vReDwggOQLEIcACKRYADUCxycEtQtoh2+WS97Ly19YMsX/q882tlU8mCYEm6cqq+KHjNynqujxQcBoURHIBiEeAAFIsAB6BY5OAgKV+Ldu7q+o/HKy6t5+BeOFU/BFOSNp7b7nBLNtZjUBjBASgLine truncated
|
||||||
|
"text/plain": [
|
||||||
|
"<Figure size 432x288 with 1 Axes>"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"metadata": {
|
||||||
|
"needs_background": "light"
|
||||||
|
},
|
||||||
|
"output_type": "display_data"
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"plt.imshow(unnormalized_posterior, cmap='Blues', aspect=aspect, extent=extent)\n",
|
||||||
|
"plt.xlabel(r'$\\mu$')\n",
|
||||||
|
"plt.ylabel(r'$\\sigma$')"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 57,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"import pymc3 as pm"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 58,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"name": "stderr",
|
||||||
|
"output_type": "stream",
|
||||||
|
"text": [
|
||||||
|
"Auto-assigning NUTS sampler...\n",
|
||||||
|
"Initializing NUTS using jitter+adapt_diag...\n",
|
||||||
|
"Sequential sampling (1 chains in 1 job)\n",
|
||||||
|
"NUTS: [sigma, mu]\n",
|
||||||
|
"100%|██████████| 10500/10500 [00:07<00:00, 1402.54it/s]\n",
|
||||||
|
"The acceptance probability does not match the target. It is 0.8838815681167924, but should be close to 0.8. Try to increase the number of tuning steps.\n",
|
||||||
|
"Only one chain was sampled, this makes it impossible to run some convergence checks\n"
|
||||||
|
]
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"with pm.Model():\n",
|
||||||
|
" # priors\n",
|
||||||
|
" mu = pm.Normal('mu', mu=200, sd=15)\n",
|
||||||
|
" sigma = pm.HalfCauchy('sigma', 10)\n",
|
||||||
|
" \n",
|
||||||
|
" # likelihood\n",
|
||||||
|
" observed = pm.Normal('observed', mu=mu, sd=sigma, observed=lengths)\n",
|
||||||
|
" # sample\n",
|
||||||
|
" trace = pm.sample(draws=10000, chains=1)"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 59,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"name": "stderr",
|
||||||
|
"output_type": "stream",
|
||||||
|
"text": [
|
||||||
|
"/opt/miniconda3/lib/python3.7/site-packages/matplotlib/axes/_base.py:3604: MatplotlibDeprecationWarning: \n",
|
||||||
|
"The `ymin` argument was deprecated in Matplotlib 3.0 and will be removed in 3.2. Use `bottom` instead.\n",
|
||||||
|
" alternative='`bottom`', obj_type='argument')\n"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAEYCAYAAABBWFftAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd4XGeV+PHvGVWrWbZluchy5BqXxLETpzghgTiE9GQ3JEtYlhAWyAYIS1sCAX4sC8vSW4AlhISeAiQshPTeEyd24l7lLluyiq3epfP7496RR6OZ0Z3RNEnn8zzzaObOvTNn7oxm7rnv+55XVBVjjDHGGGOMMSPnS3UAxhhjjDHGGDNWWIJljDHGGGOMMXFiCZYxxhhjjDHGxIklWMYYY4wxxhgTJ5ZgGWOMMcYYY0ycWIJljDHGGGOMMXFiCZYxxhhjjDHGxIklWMYYY4wxxhgTJ5ZgGWOMMcYYY0ycWIJlTJoQkX0i8jkR2SgibSJyl4hME5FHRaRFRJ4SkUki8g4RqQqx7TtTFbsxxpixzX6jjPHOEixj0su7gQuBhcAVwKPAF4ESnP/Xf09daMYYY8Y5+40yxgNLsIxJLz9R1SOqegh4EVijqm+pahfwf8CK1IZnjDFmHLPfKGM8sATLmPRyJOB6R4jbBckNxxhjjBlgv1HGeGAJljGjTxuQ578hIhnA1NSFY4wxxgyw3ygz7lmCZczosxPIFZHLRCQL+DKQk+KYjDHGGLDfKGMswTJmtFHVJuBjwJ3AIZyzhVURNzLGGGOSwH6jjAFR1VTHYIwxxhhjjDFjgrVgGWOMMcYYY0ycWIJljDHGGGOMMXFiCZYxxhhjjDHGxIklWMYYY4wxxhgTJ5mpDiCeSkpKtKKiItVhGGOMidG6devqVXVMzpljv1HGGDO6ef2NGlMJVkVFBWvXrk11GMYYY2IkIvtTHUOi2G+UMcaMbl5/o6yLoDHGGGOMMcbEiSVYxhhjjDHGGBMnY6qLoDHxsuFgI39ed5BdR1rp7OljXmkB71w8jYuWTifDJ6kOzxhjTAz6+pWHNh7m5LKJzJ1akOpwjDFjlCVYxgSob+3ilvs38sz2WvKzM1g0o4jC3Cxe2FnPX948xNySfL59zTJOr5ic6lCNMcZEqbu3H4Bdta2WYBkzCrV391LX0sUJU/JTHUpElmAZ49pU1cSNv1/L0bZubrn4RN5/1gkU5mYBzlnPx7fU8O3HtnPdHa9x6yWL+PC5c1McsTHGmGiIdUAwZlR7pbKBtu5eyoonkJmRviOdLMEyBthR08L77nyNwtwsHvjo2ZxUNnHQ/Rk+4dKTZ/C2BSXc8ueN/PfD2+jq7efj589PUcTGGGNipZrqCIwxsehyW6HT/V84fVM/Y5LkcGMH1/9qDROyM7jvxrOGJFeBinKz+Nn7TuUfls/ku4/v4Dcv701ipMYYY4wxJt1ZC5YZ13r7+vnEvW/R1tXH/R9dRfnkvGG3yfAJ37v2FNq6+/jaQ1tZOK2Qs+eXJCFaY4wx8ZHu57+NGV/21LXS0NY9Zsa4WwuWGddue6aSdfuP8Y1/PIlF04s8b5eZ4eOH71nOnJJ8br73LaqbOhIYpTHGmHiwMVjGpKdNh5o43Dh2jqUswTLj1qaqJn76zC6uPrWMq5aXRb19QU4md1y/ks6ePm65fyNqnfqNMWZUsK9rk056+/o5eLQ91WGMKul+rsQSLDMu9fUrX/rrJibn5/CfVyyN+XHmTS3gi5cu5sVd9dy95kAcIzTGGBNvkvaHZaFVHWvnb+sPDZSZN2PLhqom3jxwjKNt3akOxcRJQhMsEblYRHaISKWIfCHE/SIit7n3bxSRUwPuKxaR+0Vku4hsE5FViYzVjC/3rNnPxqom/t/li5k4IWtEj/W+M2fztvkl/M8j2zg0hpq3jTFmrBptDVi769oAZw4gM/Z09fQBzslfMzYkLMESkQzgZ8AlwBLgvSKyJGi1S4AF7uVG4OcB9/0YeExVFwGnANsSFasZX5o7e/j+kzs5e94Urjxl5ogfT0T45tUn06/K1/++NQ4RGmOSSUTKReRZ92TeFhH5pLv8u+5Jvo0i8n8iUuwurxCRDhFZ715uT+0rMGOdv93Njr/jo69f6emz1kCTOIlswToDqFTVParaDdwHXBW0zlXA79TxGlAsIjNEpAg4D7gLQFW7VbUxgbGaceQXz++msb2HL166GInTiOfyyXl8YvUCHttSw7Pba+PymMaY6InIBBE5McrNeoHPqupi4Czg4+4JwSeBk1R1GbATuDVgm92quty93BSX4I0xSfHs9loe2VSd6jBMDHSUtD8nMsEqAw4G3K5yl3lZZy5QB/xaRN4SkTtFJD/Uk4jIjSKyVkTW1tXVxS96MyYdae7krpf2ctXymRHnu4rFR86dy9yp+Xztoa12ZsyYFBCRK4D1wGPu7eUi8uBw26lqtaq+6V5vwekxUaaqT6iqv0/Wa8CsxERuks2KXIxvbdbVMqV6+vppaO0a0WOk+79wIhOsUE0Dwfsj3DqZwKnAz1V1BdAGDBnDBaCqd6jqSlVdOXXq1JHEa8aBnz1bSV+/8tkLoz3BPbzsTB9funQxe+vbuPd1K3hhTAp8Faf3RCOAqq4HKqJ5ABGpAFYAa4Lu+lfg0YDbc9wTgM+LyLkRHs9OAhozBtSPMCEwx72x9ygvVdbTO4ZPRicywaoCygNuzwIOe1ynCqhSVf8P3P04CZcxMatp6uS+1w9yzWmzmD1l+AmFY7F6USlnzZ3Mj5/aRUtnT0KewxgTVq+qNsW6sYgUAA8An1LV5oDlX8LpRni3u6gamO2eAPwMcI/btX0IOwmYnkZLNyO/8T5/V29fP03tyflNDdUD5eDRdl6urLdS6nHS1OG8l2N5TGEiE6w3gAUiMkdEsoHrgOCuGg8C17vVBM8CmtyuGjXAwYB+9BcAVj3AjMgvXthNnyofe8f8hD2HiHDrJYtpaOvmjhf2JOx5jDEhbRaRfwYyRGSBiPwEeMXLhiKShZNc3a2qfwlY/gHgcuB96k52p6pdqtrgXl8H7AYWxvelmPGgs6ePZjsZN6y1+4/x3M7ahHe/r6xt5ZFN1XR09w1a3trldClsD1puYjOG86oBCUuw3H7rNwOP4/Rn/5OqbhGRm0TEPyD4EWAPUAn8EvhYwEN8ArhbRDYCy4H/SVSsZuyrb+3injUH+McVZZRPTkzrld8p5cVcccpMfvniHo40dyb0uYwxg3wCWAp0AfcCzcCnhttInGo3dwHbVPUHAcsvBj4PXKmq7QHLp7qVchGRuTiVcO2MyjBUNW0mZE+TMHhq2xErjORBY7szP1R/HN64SKXQq5ucqVasHP7otK++jS2HY+7EEFeZiXxwVX0EJ4kKXHZ7wHUFPh5m2/XAykTGZ8aPe9ccoKu3n5vePi8pz3fLRSfy2OZqfvjkTr717mVJeU5jxjs3CfqSe4nGOcD7gU0ist5d9kXgNiAHeNKtOPqaWzHwPOBrItIL9AE3qerROLyEtNLS2cOxtp64dal+YusRVJWLT5oRl8cbCxI571Fjeze7altZecKkuFXMHQuqjh3v5tfU3kNPfz8lBTmD1ukLk8h52Y3P7ails6efi0+aPqI4x7J4fBpDvUUbqpyC40tnxreIWSwSmmAZkw56+vr5/Wv7OW/hVOaXFiTlOcsn53H9qgp+/fJePnzuHOaXFibleY0Zz0TkWUL0PlHV1ZG2U9WXCP2b/0iIZajqAzjdCce053fW0devcUuwOnvSp3tVmjRgJdQb+47R3t1L+4wi8nPGxuFevFsen9vptB5etdwpci3u18Cruxu49OQZZGVE39HLP74oGuPh8zjeJHIMljFp4dHNNdS2dPHBsyuS+rw3nz+fvOxMfvDkzqQ+rzHj2H8An3Mv/w+nZPvalEY0iiWydcWYdNPU0TOoa2BXb/Iq3Pn/18bL/1y8X2U67jdPCZaInJToQIxJlN+8vJeKKXm8fWFyK3hNys/mQ2+bwyObath8KD36BBszlqnquoDLy6r6GeDMVMc1Hjy97Qiv7m5IdRhmzEp8F8fndtTSkaJW1mPuGLNdtS0pef7RJDiZ6uju46GNh9lb3xZy/Tf2HeVAQ/KrP3ptwbpdRF4XkY+JSHFCIzImjjYcbOTNA4184OwKfL7k90H/8LlzKM7L4ntP7Ej6cxsz3ojI5IBLiYhcBNhAiCRo7eqltmVsF/Vp7uzhb+sPUZuC4kX+Y0orvjBy6VLgJJ4e2VTNvjAJxljV2OEkpf4Kj4cbO0Kud7ixg7cOHktaXH6eEixVfRvwPpw5q9aKyD0icmFCIzMmDn77yj7yszO45rRZKXn+wtwsbnr7PJ7bUcfafWNuDLwx6WYdTpfAdcCrwGeBD6U0onFuw8FG/rb+UKrDCKmrd2hrRXVTB39bf4jK2lb+tv7QoHLdx9qcA7rDTclPsPxV9NYfbEz6c49UU3sP/WnYhStVjrZ18+yO2pDd2iSGljpVpaevn42jqKdMPE53t3X18crueo66/5fpljh7HoOlqruAL+OUrH07cJuIbBeRqxMVnDEjUdvSyd83HubaleUU5malLI4PrKpgamEO3318R9qUKDZmLFLVOao61/27QFXf5RawMAlQ29I5bIvKvob0PKuuqjy2uWbIVBoHjzpnwf2lngPnqGrrcpKt/Q1t7E/T15VuWrt6eW5nLVsONw+/chob7qdbVak61u7pN37L4SaaO3pCFsPIyYq9NEI0xxfxOBbp7evnqa1HqGnqDNt6FPb5R/zsTtGcupYuttc0u48Z+VGTffzldQzWMhH5Ic58VquBK1R1sXv9hwmMz5iY3bvmID19yvWrTkhpHBOyM7j5/Pms2XuUlyrrUxqLMWORiFwd6ZLq+MaqV3c38Mwon8PptT2Rx40FnmkPHB+TyJakju4+XtxVR3cSiiw0dfTwxr6jcWthqjrWzp661oHLine truncated
|
||||||
|
"text/plain": [
|
||||||
|
"<Figure size 864x288 with 4 Axes>"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"metadata": {
|
||||||
|
"needs_background": "light"
|
||||||
|
},
|
||||||
|
"output_type": "display_data"
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"pm.traceplot(trace);"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 60,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"text/plain": [
|
||||||
|
"(0, 15.0)"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"execution_count": 60,
|
||||||
|
"metadata": {},
|
||||||
|
"output_type": "execute_result"
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7oAAAF+CAYAAACs4mbmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3X2QZNlZ3/nvc869mdXdMxoJzQDSzEhgW4YwAgQMwjY2KxtkizFCgdcblvArmBjjRcGySxBgvGFkw653wy/hwLKRBcgYAxJ+k1c2I17sheXFvGikFUISb4OQrNFoJQ2geemuyrz3nGf/OOdmZddUVddMn6rKzv59Qq3pqrx58ubNl77PfZ7zHHN3RERERERERLZFOO8dEBEREREREWlJga6IiIiIiIhsFQW6IiIiIiIislUU6IqIiIiIiMhWUaArIiIiIiIiW0WBroiIiIiIiGwVBboiIiJyTWb2EjN7z3nvh4iIyEko0BUREVljZk+s/clmtrv281847/17Oszsa8zsp65nDHf/KXf/jEa7JCIicqq6894BERGRTeLut0x/N7P3A1/j7v/5qO3NrHP38Sz27bxc73M0s+juqeU+iYiIHEcZXRERkafAzL7DzH7YzN5oZo8Df9HMfsDMXrO2zZfUIHn6+S4ze7OZfczMftvMvu6Y8X/AzP6pmf0XM3vczH7SzO5eu/2PmdkDZvaomf2SmX3B2m1/zczeX+/3PjN7pZl9JvBa4I/XrPQjddsdM/tHZvZBM/uImf0zM9tZ338z+1Yz+/+A7z7kOX2Gmf0/ZvZxM/sVM/szhzyHHzWzy8Afv55jLiIi8lQp0BUREXnqvgL4IeA24IeP29DMIvCfgLcBdwIvBb7JzL74mLv9ReBvA7cD7wX+VR3rduBHgH8IPBv4TuB+M3uWmT0D+EfAS939VuALgXe5+68ArwZ+xt1vcffb62P8A+BTgc8CXgB8CvC31vbhLuAW4HnA/3jgOc3qc/oR4A7gfwZ+2Mz+wNpmXwn8HeBW4OePO0YiIiKtKdAVERF56n7W3f+ju2d3373Gtn8YeIa7/+/uvnT3B4HvBV55zH3+o7v/nLsvgG8FvsjMngO8HHiPu7/R3Ud3/wHgfcCUTXXghWa24+4fdvf3Hja4mQXga4BvcPffc/fHgL93YJ9G4DV1nw8+xy8EZsDfd/ehlna/9cD93+zuP1+P0eIax0hERKQpzdEVERF56j74FLZ9PvA8M/v42u8i8FMnGd/dHzWzR4Hn1j8fOLDtB4A73f0xM3sV8I3AvzCznwX+F3f/jUPG/2RgDvyymU2/swPbfMTdl0fs33OB/+bufnA/DnsOIiIiZ00ZXRERkafOD/x8Gbi49vMnr/39g8Bvuvsz1/7c6u4vP2b89Tm5t1FKpB+uf55/YNvnAR8CcPe3uvuXAM8BHgT++RH7+xFgCXza2j7d5u63HfMc1z0M3G1rUfL6fpzg/iIiIqdKga6IiMj1eyfwZ+pc2ecAX792288DSzP7xtoAKprZZ5rZ5x0z3svN7I+Y2Rz4Dkqp9Icp82I/w8z+vJl1ZvaVwB+gzNN9jpm93MwuUoLYy8DU6fgjwF1m1gPUDsjfA/xjM7vDirvM7E+d8Pn+V0pp8zeaWW9mfxK4F/jXJ7y/iIjIqVKgKyIicv2+D/hVSvnujwJvmm6oy/LcC7wYeD/wCCXT+oxjxvsBSoD7CKVZ1F+qY30M+HLgm4HfoTSB+jJ3/11KOfQ3AR+ut/1RShMqgJ8AfhP4SO2iDKXE+QPALwGPAj9OaUp1TXXO7cuBV9R9/E7gK48okxYRETlzdvX0GhERETlPZvYDwIPu/prz3hcREZEblTK6IiIiIiIislVOLdA1szeY2UfN7N1rv3uNmX3IzN5Z/9x7xH1fZma/bmYPmtm3nNY+ioiIiIiIyPY5tdJlM/si4Ang+939hfV3rwGecPd/cMz9IvAbwEuBh4C3Aa86ai1AERERERERkXWnltF1958Gfvdp3PXFlLlJ76vr972J0uxCRERERERE5JrOY47uq83sXbW0+VmH3H4nVy8y/xBXL0AvIiIiIiIicqTujB/vu4Bvpywi/+3APwS++sA2dvBOHLPovJndB9wHcOnSpc/7tE/79DZ7KnKDUz/147WateEbdKQ3rYn+pu1PM4f9K3V+wzRhzZ7TZh2cTTrGm0bHRkRuRB/4wPt55JFHTvQVdqaBrrt/ZPq7mX03ZeH7gx4C7l77+S7g4WPGfD3weoDP+7x7/Od+8YE2Oytyg2s1/75VsJI3bH9So4FS3pxxWu1LbjVOq4sJrV70RtFcaBQhhAb7ExvtTNdonNBonFbPq9VrZRv23mmh1XMSETlLX/gF95x42zMtXTaz56z9+BXAuw/Z7G3AC8zsU81sBrwSeMtZ7J+IiIiIiIjc+E4to2tmbwReAtxuZg8B3wa8xMxeRKmqfD/w1+u2zwW+x93vdffRzF4N/BgQgTe4+3tOaz9F5Ca1YWW1LZKW7bLvbcZpVlXQZBSaHaDcqOhzkxJqG/ZxEBERuW6nFui6+6sO+fX3HrHtw8C9az/fD9x/SrsmIiIiIiIiW+ysm1GJiFyXjcs8bdAOtdqVVuM0ywy3GaaZ2GiPNimLv2mtiTZtXraIiNx4zmN5IREREREREZFTo0BXREREREREtopKl0XkprRpZb4tSjXbLSnVanmhzSnxhXZVrM2aUbUoXb7+Ico4zWqgVSp8FrQ0kIjItSmjKyIiIiIiIltFga6IiIiIiIhsFZUui4hch01aK7ZVl+NNG6fVMbZmZbVt9ic02J92XZfb2NZuyZu1NyIichLK6IqIiIiIiMhWUUZXRG5Km5SJhUbZz1ZNpBqlYts1x2oyDK1eLW+U32vSgKzVvjQZRUREZHMooysiIiIiIiJbRYGuiIiIiIiIbBWVLovIDaXdep9tNNudBgO1ayLVZqDUaodaabU7jS4RuzcoO97S9W+1TqyIiFwvZXRFRERERERkqyjQFRERERERka1yaqXLZvYG4MuAj7r7C+vv/j7wcmAJ/BbwVe7+8UPu+37gcSABo7vfc1r7KSI3pw0rqm1Sdtyq5Hhb19ENjcphW1ULb9J7cMNmBIiIiFy308zofh/wsgO/+wnghe7+WcBvAH/zmPv/CXd/kYJcEREREREReSpOLaPr7j9tZp9y4Hc/vvbjLwB/7rQeX0TkLGzSeryblkFtth5vk1HajWQbdJxbraMrIiKybc5zju5XA2894jYHftzM3m5m9x03iJndZ2YPmNkDH3vkY813UkRERERERG4s5xLomtnfAkbgB4/Y5Avd/XOBLwW+zsy+6Kix3P317n6Pu99zx+13nMLeioiIiIiIyI3kzNfRNbO/QmlS9cV+RN2Wuz9c//tRM3sz8GLgp89uL0VkUzUrY92whkJNylg3bP3bVs2xmjV/alTl26qpVYvDrCZSIiIihzvTjK6ZvQz4ZuDL3f3KEdtcMrNbp78Dfwp499ntpYiIiIiIiNzITi3QNbM3Aj8PfJqZPWRmfw14LXAr8BNm9k4ze13d9rlmdn+96ycBP2tmvwz8EvAj7v6jp7WfIiIiIiIisl1Os+vyqw759fcese3DwL317+8DPvu09ktEBFqWHLcZp0UZa6uS400bp9UxDqFVyXGr59Wi67LcSKxR2buIiFzbeXZdFhEREREREWnuzJtRiYhsk01qRtVqHd1W42xaRreV3Cgz3OR5bdrB2TCtMqhKxIqI3HiU0RUREREREZGtokBXREREREREtopKl0XkptRqzdlW47So8s0b1kRqW0uXY6vn1WSUzdKsVLjJKCIicjNTRldERERERES2igJdERERERER2SoqXRYRuQ6tymrbdF1uszNjyk3GGVKjsu5GRb6zRtd2tY6uiIjI5lNGV0RERERERLaKAl0RERERERHZKipdFpEbS6tS4TbDNBunRYfiTeuWvMipyTitysNDo16+fbNy9c0YYxM1at6s7s0iIjcxZXRFRERERERkqyijKyI3pU1qIlXGuf4xWmVil7lNM6rF2GacVmKjNOE8b1ZTq02iDKqIiGyKU8vomtkbzOyjZvbutd99gpn9hJn9Zv3vs46478vM7NfN7EEz+5bT2kcRERERERHZPqdZuvx9wMsO/O5bgP/i7i8A/kv9+SpmFoF/Cnwp8IeAV5nZHzrF/RQREREREZEtcmqly+7+02b2KQd+/QrgJfXv/xL4KeCbD2zzYuBBd38fgJm9qd7vvae0qyJbyRqVabYqzd1WrQ5PizLWVqXLe41KjvdSm2ZUrYRGn4mL+kzcdFp9n4qIyNk562ZUn+TuHwao//3EQ7a5E/jg2s8P1d8dyszuM7MHzOyBjz3ysaY7KyIiIiIiIjeeTey6fNhl0yMvn7v76939Hne/547b7zjF3RIREREREZEbwVl3Xf6ImT3H3T9sZs8BPnrINg8Bd6/9fBfw8JnsnYjIU7RJ6+gOqdH6t41KjluN00oX2pSfptzmn86trIDesApfVRyLiNy8zjqj+xbgr9S//xXg/zpkm7cBLzCzTzWzGfDKej8RERERERGRazq1jK6ZvZHSeOp2M3sI+Dbg/wD+tZn9NeC/Af9D3fa5wPe4+73uPprZq4EfAyLwBnd/z2ntp4gcr1lGpFX2asMyNK2adbXoIzU0Wv+2VROpy0ObcRolYulDm2u7Y7PXfHNSuq0+5xv28RQRkZvYaXZdftURN33xIds+DNy79vP9wP2ntGsiIiIiIiKyxTaxGZWIiIiIiIjI03bWzahERLZKs3V0G9QuL1Kj0uWxUenyss3+tCpd3una7E+r9Yo3yaaVHGvdWhERuV7K6IqIiIiIiMhWUaArIiIiIiIiW0WlyyIi16FVEWuLDryt1q19fNlmnMf22ozTxzZlrLfM25Quj95mnBaadUtuNFCzcZqMIiIiNzNldEVERERLine truncated
|
||||||
|
"text/plain": [
|
||||||
|
"<Figure size 1152x432 with 2 Axes>"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"metadata": {
|
||||||
|
"needs_background": "light"
|
||||||
|
},
|
||||||
|
"output_type": "display_data"
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"fig, axes = plt.subplots(2, sharex=True, sharey=True, figsize = (16, 6))\n",
|
||||||
|
"axes[0].imshow(posterior, cmap='Blues', extent=extent, aspect=1)\n",
|
||||||
|
"axes[0].set_ylabel('$\\sigma$')\n",
|
||||||
|
"axes[1].scatter(trace['mu'], trace['sigma'], alpha=0.01)\n",
|
||||||
|
"axes[1].set_ylabel('$\\sigma$')\n",
|
||||||
|
"axes[0].set_title('True posterior')\n",
|
||||||
|
"axes[1].set_title('Sampled $\\\\theta$')\n",
|
||||||
|
"plt.xlabel('$\\mu$')\n",
|
||||||
|
"plt.xlim(150, mm.max())\n",
|
||||||
|
"plt.ylim(0, ss.max())"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 4,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"import pyro\n",
|
||||||
|
"import pyro.distributions as dist\n",
|
||||||
|
"import pyro.optim\n",
|
||||||
|
"from pyro.infer import SVI, Trace_ELBO\n",
|
||||||
|
"import torch\n",
|
||||||
|
"import torch.distributions.constraints as constraints"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 38,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"\n",
|
||||||
|
"def model():\n",
|
||||||
|
" # priors\n",
|
||||||
|
" mu = pyro.sample('mu', dist.Normal(loc=torch.tensor(200.), \n",
|
||||||
|
" scale=torch.tensor(15.)))\n",
|
||||||
|
" sigma = pyro.sample('sigma', dist.HalfCauchy(scale=torch.tensor(10.)))\n",
|
||||||
|
" \n",
|
||||||
|
" # likelihood\n",
|
||||||
|
" with pyro.plate('plate', size=2):\n",
|
||||||
|
" pyro.sample(f'obs', dist.Normal(loc=mu, scale=sigma), \n",
|
||||||
|
" obs=torch.tensor([195., 185.]))\n",
|
||||||
|
" \n",
|
||||||
|
"def guide():\n",
|
||||||
|
" # variational parameters\n",
|
||||||
|
" var_mu = pyro.param('var_mu', torch.tensor(180.))\n",
|
||||||
|
" var_mu_sig = pyro.param('var_mu_sig', torch.tensor(5.), \n",
|
||||||
|
" constraint=constraints.positive)\n",
|
||||||
|
" var_sig = pyro.param('var_sig', torch.tensor(5.))\n",
|
||||||
|
" \n",
|
||||||
|
" # factorized distribution\n",
|
||||||
|
" pyro.sample('mu', dist.Normal(loc=var_mu, scale=var_mu_sig))\n",
|
||||||
|
" pyro.sample('sigma', dist.Chi2(var_sig))\n",
|
||||||
|
" "
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 53,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"name": "stdout",
|
||||||
|
"output_type": "stream",
|
||||||
|
"text": [
|
||||||
|
"[iteration 1] loss: 20.5714\n",
|
||||||
|
"[iteration 101] loss: 12.0774\n",
|
||||||
|
"[iteration 201] loss: 11.7402\n",
|
||||||
|
"[iteration 301] loss: 9.8967\n",
|
||||||
|
"[iteration 401] loss: 28.1990\n",
|
||||||
|
"[iteration 501] loss: 13.2090\n",
|
||||||
|
"[iteration 601] loss: 11.4618\n",
|
||||||
|
"[iteration 701] loss: 14.9521\n",
|
||||||
|
"[iteration 801] loss: 11.2410\n",
|
||||||
|
"[iteration 901] loss: 9.7651\n",
|
||||||
|
"[iteration 1001] loss: 9.7140\n",
|
||||||
|
"[iteration 1101] loss: 9.5427\n",
|
||||||
|
"[iteration 1201] loss: 9.0695\n",
|
||||||
|
"[iteration 1301] loss: 9.0411\n",
|
||||||
|
"[iteration 1401] loss: 10.2914\n",
|
||||||
|
"[iteration 1501] loss: 8.6590\n",
|
||||||
|
"[iteration 1601] loss: 9.7012\n",
|
||||||
|
"[iteration 1701] loss: 12.4810\n",
|
||||||
|
"[iteration 1801] loss: 8.0311\n",
|
||||||
|
"[iteration 1901] loss: 9.6596\n",
|
||||||
|
"[iteration 2001] loss: 10.1730\n",
|
||||||
|
"[iteration 2101] loss: 8.9640\n",
|
||||||
|
"[iteration 2201] loss: 9.5495\n",
|
||||||
|
"[iteration 2301] loss: 9.6852\n",
|
||||||
|
"[iteration 2401] loss: 10.8672\n",
|
||||||
|
"[iteration 2501] loss: 8.9980\n",
|
||||||
|
"[iteration 2601] loss: 8.3449\n",
|
||||||
|
"[iteration 2701] loss: 9.4454\n",
|
||||||
|
"[iteration 2801] loss: 9.4535\n",
|
||||||
|
"[iteration 2901] loss: 9.0612\n",
|
||||||
|
"[iteration 3001] loss: 8.6264\n",
|
||||||
|
"[iteration 3101] loss: 8.5423\n",
|
||||||
|
"[iteration 3201] loss: 7.1836\n",
|
||||||
|
"[iteration 3301] loss: 8.9224\n",
|
||||||
|
"[iteration 3401] loss: 8.8739\n",
|
||||||
|
"[iteration 3501] loss: 9.0185\n",
|
||||||
|
"[iteration 3601] loss: 38.2326\n",
|
||||||
|
"[iteration 3701] loss: 8.8873\n",
|
||||||
|
"[iteration 3801] loss: 9.0606\n",
|
||||||
|
"[iteration 3901] loss: 8.0937\n",
|
||||||
|
"[iteration 4001] loss: 9.1192\n",
|
||||||
|
"[iteration 4101] loss: 7.1300\n",
|
||||||
|
"[iteration 4201] loss: 9.4013\n",
|
||||||
|
"[iteration 4301] loss: 8.7548\n",
|
||||||
|
"[iteration 4401] loss: 9.2406\n",
|
||||||
|
"[iteration 4501] loss: 9.3022\n",
|
||||||
|
"[iteration 4601] loss: 8.9886\n",
|
||||||
|
"[iteration 4701] loss: 9.1033\n",
|
||||||
|
"[iteration 4801] loss: 9.0818\n",
|
||||||
|
"[iteration 4901] loss: 8.5367\n"
|
||||||
|
]
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"pyro.clear_param_store()\n",
|
||||||
|
"pyro.enable_validation(True)\n",
|
||||||
|
"\n",
|
||||||
|
"svi = SVI(model, guide, \n",
|
||||||
|
" optim=pyro.optim.ClippedAdam({\"lr\":0.01}), \n",
|
||||||
|
" loss=Trace_ELBO())\n",
|
||||||
|
"\n",
|
||||||
|
"# do gradient steps\n",
|
||||||
|
"c = 0\n",
|
||||||
|
"for step in range(5000):\n",
|
||||||
|
" c += 1\n",
|
||||||
|
" loss = svi.step()\n",
|
||||||
|
" if step % 100 == 0:\n",
|
||||||
|
" print(\"[iteration {:>4}] loss: {:.4f}\".format(c, loss))"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 54,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"name": "stdout",
|
||||||
|
"output_type": "stream",
|
||||||
|
"text": [
|
||||||
|
"190.1922149658203 5.435718059539795\n"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"text/plain": [
|
||||||
|
"10.607619285583496"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"execution_count": 54,
|
||||||
|
"metadata": {},
|
||||||
|
"output_type": "execute_result"
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"print(pyro.param('var_mu').item(), pyro.param('var_mu_sig').item())\n",
|
||||||
|
"pyro.param('var_sig').item()"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 55,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"sigma = dist.Chi2(pyro.param('var_sig')).sample((10000,)).numpy()\n",
|
||||||
|
"mu = dist.Normal(pyro.param('var_mu'), pyro.param('var_mu_sig')).sample((10000,)).numpy()"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 56,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"text/plain": [
|
||||||
|
"(0, 15.0)"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"execution_count": 56,
|
||||||
|
"metadata": {},
|
||||||
|
"output_type": "execute_result"
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7oAAAGACAYAAACDeh6PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xm8ZHlZ5/nPc7ZY7pZrVWXWAsggKIUUUqDd9iiKKFQLjD12D2gLrTilvqBtdKYF0VeLQrf2jPY4Nq1YCDI4LNraOoWWArbSqKNAFQMICFoUBZW15nqX2M72zB/n3Kyoy72ZWZW/u2Tk9/163VdGxDnxi1+ciLh5n3ie8/zM3RERERERERGZFdFuT0BEREREREQkJAW6IiIiIiIiMlMU6IqIiIiIiMhMUaArIiIiIiIiM0WBroiIiIiIiMwUBboiIiIiIiIyUxToioiIiIiIyExRoCsiIiIiIiIzRYGuiIjIjDKzT5vZc3bosX7OzF59gfu+3sxe/xge4yNm9tRHPTkREbnsKNAVEZFLnpmtTf3UZjaauv49uz2/x8LM7jazb72YMdz9qe7+wUDz2W9m3h7ToZndtx7Ymtlh4GXAr224zzPM7C/b/T9iZtdd4PhrZvYlM/sdM7t+ardfAH42xPMREZHZpkBXREQuee4+v/4DfAl44dRt75ze18yS3ZnlzrmY53iO+94AHG+PaR/4YeD/MLNrgH8B3Obuo6lxrgFuA/49cBC4C/ipczz0DcCpqdfxGcAngA+b2VPafW4FvtnMjjzW5yciIpcHBboiIjLz2uzoa8zsk8DAzJI2e/jfTe3zdjN7Y3v5qJn9rpkdN7MvmNmPXMD4P2FmnzGz02b2G2bWbbd9lZl90MzOtKXEL9pw39eY2b1mtmpmnzOz55rZbwLXAe9ts5s/fr45bfEcz2aFzzWPze67ydO8Afjo1PUPt/9mwAuA/7Zh/18E3uLut7YB8HuAZ53jMN4AfHz9irufdPc3AB8DXtHeNgbuAL7tHOOIiIgo0BURkcvGS4F/DOxz93KrncwsAt5Lk028Gngu8Goz+/bzjP89wLcDTwS+EvgpM0vbsd4PXAH8S+CdZvbk9rGeDLwKeJa7L7T3v9vdv5epzDRNye6FzGnT53i+eVzg8XkG8JF2vH3Av6UJOr8APA343NTjLQIvBn596v4RMD7H8XsGU4HulM+2z3nd3wJPP8c4IiIiCnRFROSy8cvufs90ee0WngUcdvefdffc3e8C3gK85Dz3e1M7/imaIPClwNcD88DPt2P9KfAH7TaACugAX21mqbvf7e6fv4g5bfUczzePc9133Q3AvzazUzQBrtME4g7sA1an9n0ukAKfbDPIZ4B3Al/cYuz18TcLdJeA41PXV9vHExER2dLMn6ckIiLSuucC93sccLQNztbFwJ8/ivG/CBxtf+5x93rDtqsB3P3OtqHT64Gnmtn7gB9z9/se45y2eo7nnMd57ouZdYCvAp7g7sc22eU0sDB1/fHAre7+nVNj3Ab88XnG/8SG22PgHwLTZdoLwPRxEBER+TLK6IqIyOXCN1wfAv2p61e1/94DfMHd9039LLj7TecZ/9qpy9cB97U/17bl0NPb7j07Kfd3ufs/oglmnaZ508b5XuicNj7HdeedxznuC3A9MNgiyAX4JE259roOzfEFwMyeANxI00xqq/FrmrLkaT8E5DRl1+u+LCAWERHZSIGuiIhcrj4OfLeZxWb2fOCb2ts/Aqy0zZl67fbrzexcjZQAXmlm15jZAeB1wG/RNGwaAD9uZqk1a9q+kKYxE2b2ZDP7ljajOQZGNOXMAA8CX3GRc1p3znlcgGcAnz7H9tt4+PhB07Tqm9oGWtcC7wJ+si3r3nJ8dy8AzOzatjHYzwAvmbq9AzwT+MAFzltERC5TCnRFRORy9a9ogr0zNI2kfh/A3av29htoGi2doGmqtHSe8d5F0+zprvbnje6eAy+i6Up8AvgV4GXu/tn2Ph3g59ttD9A0inpdu+3naBpanQF+9DHOifY5nW8e53MD8KlzbH8HcJOZ9drrf0qThf074C+A33T3t5xn/K9pO0+fBv4E2A/c6O4fmdrvRcAHNyntFhEReQRrekiIiIjIY2VmdwM/4O5/sttz2S1m9u+Ah9z9ly5g39cDuPvrH+VjfBh4hbufK+gWERFRMyoRERG5eO7+uvPvddGP8XXb/RgiIjIbtq102czeZmYPmdmnpm57vZnda2Yfb382bexhZs83s8+Z2Z1m9trtmqOIiIjsig+2PyIiItti20qXzewbgTXgHe5+fXvb64E1d/+Fc9wvpjmn53nAMZqGFi91989sy0RFRERERERkpmxbRtfdPwRs1V3xXJ4N3Onud7XNM94DvDjo5ERERERERGRm7UbX5VeZ2Sfb0ub9m2y/mkcuWn+MRy5oLyIiIiIiIrKlnW5G9avAG2gWpX8D8IvA92/Yxza535b11WZ2M3AzwNzc3DOf/OSnhJmpyCVO/dTPLdRZG76HjvRea6K/1+YTzGb/S+3eMEFYsOe0tw7OXjrGe42OjYhcir74xbs5ceLEBf0K29FA190fXL9sZm8B/mCT3Y4B105dvwbYcr08d78FuAXgmc+80f/yw7eHmazIJS7U+fehgpV6j82nCjRQVe+dcULNpQ41TqgvE0K96IGiuShQhBAFmE8caDJJoHGiQOOEel6hXivbY++dEEI9JxGRnfQNX3fjBe+7o6XLZnZk6up3svni8x8FnmRmTzCzDHgJcOtOzE9EREREREQufduW0TWzdwPPAQ6Z2THgp4HnmNkNNFWVdwM/2O57FPh1d7/J3UszexXwPiAG3ubun96ueYrIZWqPldWGSFqGy76HGSdYVUGQUQh2gOpARZ97KaG2xz4OIiIiF23bAl13f+kmN791i33vA26aun4bcNs2TU1ERERERERm2E43oxIRuSh7LvO0hyYUaiqhxgmWGQ4zTDBxoBntpSz+XmtNtNfOyxYRkUvPbiwvJCIiIiIiIrJtFOiKiIiIiIjITFHpsohclvZamW+IUs1wS0qFWl5o75T4Qrgq1mDNqEKULl/8EM04wWqgVSq8E7Q0kIjI+SmjKyIiIiIiIjNFga6IiIiIiIjMFJUui4hchL20VmyoLsd7bZxQx9iCldWGmU8UYD7hui6HMavdkvfWbERE5EIooysiIiIiIiIzRRldEbks7aVMLATKfoZqIhUoFRuuOVaQYQj1anmg/F6QBmSh5hJkFBERkb1DGV0RERERERGZKQp0RUREREREZKaodFlELinh1vsMI9h0AgwUrolUmIGqUBMKJdR0An1F7B6g7HhG17/VOrEiInKxlNEVERERERGRmaJAV0RERERERGbKtpUum9nbgO8AHnL369vb/nfghUAOfB74Pnc/s8l97wZWgQoo3f3G7ZqniFye9lhRbZCy41Alx7O6jm4UqBw2VLXwXnoP7rEzAkRERC7admZ03w48f8NtHwCud/evAf4O+Ilz3P+b3f0GBbkiIiIiIiLyaGxbRtfdP2Rmj99w2/unrv418F3b9fgiIjthL63Hu9cyqMHW4w0ySriRbA8d51Dr6IqIiMya3TxH9/uBP9pimwPvN7M7zOzmcw1iZjeb2e1mdvvxE8eDT1JEREREREQuLbsS6JrZTwIl8M4tdvkGd/9a4AXAK83sG7cay91vcfcb3f3Gw4cOb8NsRURERERE5FKy4+vomtnLaZpUPde3qNty9/vafx8ys98Dng18aOdmKSJ7VbAy1j3WUChIGeseW/82VHOsYM2fAlX5hmpqFeIwq4mUiIjI5nY0o2tmzwdeA7zI3Ydb7DNnZgvrl4FvAz61c7MUERERERGRS9m2Bbpm9m7gr4Anm9kxM3sF8CZgAfiAmX3czN7c7nvUzG5r73ol8Bdm9gngI8Afuvsfb9c8RUREREREZLZsZ9fll25y81u32Pc+4Kb28l3A07drXiIiELLkOMw4IcpYQ5Uc77VxQh3jKApVchzqeYXouiyXEgtU9i4iIue3m12XRURERERERILb8WZUIiKzZC81owq1jm6ocfZaRjeUOlBmOMjz2msHZ48JlUFVIlZE5NKjjK6IiIiIiIjMFAW6IiIiIiIiMlNUuiwil6VQa86GGidElW+9x5pIzWrpchzqeQUZZW8JViocZBQREbmcKaMrIiIiIiIiM0WBroiIiIiIiMwUlS6LiFyEUGW1Ybouh5lMWdVBximqQGXdgYp8s0Df7WodXRERkb1PGV0RERERERGZKQp0RUREREREZKaodFlELi2hSoXDDBNsnBAdivdat+RJXQUZJ1R5eBSol28arFx9b4yxFwVq3qzuzSIilzFldEVERERERGSmKKMrIpelvdREqhnn4scIlYnN6zDNqCZlmHFCiQOlCTv13mpqtZcogyoiInvFtmV0zextZvaQmX1q6rYDZvYBM/v79t/9W9z3+Wb2OTO708xeu11zFBERERERkdmznaXLbweev+G21wL/1d2fBPzX9vojmFkM/CfgBcBXAy81s6/exnmKiIiIiIjIDNm20mV3/5CZPX7DzS8GntNe/r+ADwKv2bDPs4E73f0uADN7T3u/z2zTVEVmkgUq0wxVmjurQh2eEGWsoUqXx4FKjsdVmGZUoUSBPhN9fSYuO6F+n4qIyM7Z6WZUV7r7/QDtv1dsss/VwD1T14+1t23KzG42s9vN7PbjJ44HnayIiIiIiIhcevZi1+XNvjbd8utzd7/F3W909xsPHzq8jdMSERERERGRS8FOd11+0MyOuPv9ZnYEeGiTfY4B105dvwa4b0dmJyLyKO2ldXSLKtD6t4FKjkONE0oShSk/reow/3XOZAX0HqvwVcWxiMjla6czurcCL28vvxz4fzbZ56PAk8zsCWaWAS9p7yciIiIiIiJyXtuW0TWzd9M0njpkZseAnwZ+HvhtM3sLine truncated
|
||||||
|
"text/plain": [
|
||||||
|
"<Figure size 1152x432 with 2 Axes>"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"metadata": {
|
||||||
|
"needs_background": "light"
|
||||||
|
},
|
||||||
|
"output_type": "display_data"
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"fig, axes = plt.subplots(2, sharex=True, sharey=True, figsize = (16, 6))\n",
|
||||||
|
"axes[0].imshow(posterior, cmap='Blues', extent=extent, aspect=1)\n",
|
||||||
|
"axes[0].set_ylabel('$\\sigma$')\n",
|
||||||
|
"axes[1].scatter(mu, sigma, alpha=0.05)\n",
|
||||||
|
"axes[1].set_ylabel('$\\sigma$')\n",
|
||||||
|
"axes[0].set_title('True posterior $P(\\\\theta|D)$')\n",
|
||||||
|
"axes[1].set_title('Approximated posterior $Q(\\\\theta)$')\n",
|
||||||
|
"plt.xlabel('$\\mu$')\n",
|
||||||
|
"plt.xlim(150, mm.max())\n",
|
||||||
|
"plt.ylim(0, ss.max())"
|
||||||
|
]
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"metadata": {
|
||||||
|
"kernelspec": {
|
||||||
|
"display_name": "Python 3",
|
||||||
|
"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.7.3"
|
||||||
|
}
|
||||||
|
},
|
||||||
|
"nbformat": 4,
|
||||||
|
"nbformat_minor": 2
|
||||||
|
}
|
||||||
Reference in new issue
Block a user