deep-learning/Models_Algorithms_More/Probabilitic Models - Bayesian Linear Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# #1 Bayesian Linear Regression\n",
    "\n",
    "Bayesian Linear Regression with toy data in 1D.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXmYHFW9//86Vb3NdM9M9hAIEFBZs5Ow3EACKgRlUVFE\nBCUoenH5uv0uileveHF58MIPkEX4IVfBFb9clesVBAyGVQ2CIhfZYjCEkJBMlpnpvbuqzu+P6qqu\nqq7qZaZ7pifU+3nmmZnuqnNObe/6nPdnOUJKSYgQIUKE2HugTPQAQoQIESJEexESe4gQIULsZQiJ\nPUSIECH2MoTEHiJEiBB7GUJiDxEiRIi9DCGxhwgRIsRehpDYQ4QIEWIvQ0jsIUKECLGXIST2ECFC\nhNjLEJmITmfMmCHnzZs3EV2HCBEixKTFk08+uVNKObPRdhNC7PPmzeOJJ56YiK5DhAgRYtJCCPFy\nM9uFUkyIECFC7GUIiT1EiBAh9jKExB4iRIgQexkmRGP3Q7lcZsuWLRQKhYkeSohxQCKRYO7cuUSj\n0YkeSogQex26hti3bNlCX18f8+bNQwgx0cMJ0UFIKdm1axdbtmzhoIMOmujhhAix16FrpJhCocD0\n6dNDUn8dQAjB9OnTw9lZiBAdQtPELoT4nhBihxDiGcdn04QQvxVCbKj8njqWwYSk/vpBeK1DhOgc\nWrHYbwNO9Xx2KfCAlPJNwAOV/0OE2CtRKOsUNX2ihxEiREM0TexSyoeB3Z6P3wHcXvn7duCdbRrX\nhEBVVRYvXsz8+fM544wzGBoamughufD2t7+9LWP66le/ylVXXQXAV77yFdauXTvmNl8P2JUtsSdb\nnuhhhAjREGPV2GdLKbcBVH7PGvuQJg49PT089dRTPPPMM0ybNo0bb7yxLe1qmtaWdu655x6mTJnS\nlrYsXH755bz1rW9ta5t7K6SUGOHi7yEmAcbNeSqE+KgQ4gkhxBODg4Pj1e2ocdxxx/Hqq6/a/195\n5ZUsX76chQsXctlll9mff+1rX+Owww7j5JNP5txzz7Ut4RNPPJF//dd/ZdWqVXz7299mcHCQd7/7\n3Sxfvpzly5fz2GOPAfDQQw+xePFiFi9ezJIlS0in02zbto2VK1fas4dHHnkEMEsx7Ny5E4Crr76a\n+fPnM3/+fK699loANm3axOGHH85HPvIRjjzySE455RTy+Xzd41yzZg3/9V//Zbd/2WWXsXTpUhYs\nWMDzzz8PQDab5UMf+hDLly9nyZIl/Pd//3c7TvGkg5QQ0nqIyYCxhjtuF0LMkVJuE0LMAXYEbSil\nvAW4BWDZsmX1n48bboC//32MQ/PgjW+ET36yqU11XeeBBx7gwx/+MAD3338/GzZs4PHHH0dKyZln\nnsnDDz9Mb28vP//5z/nLX/6CpmksXbqUo446ym5naGiIhx56CID3v//9fPazn+X4449n8+bNrF69\nmueee46rrrqKG2+8kRUrVpDJZEgkEtxyyy2sXr2aL33pS+i6Ti6Xc43vySef5Pvf/z7r169HSskx\nxxzDqlWrmDp1Khs2bOCnP/0p3/3ud3nve9/Lz3/+c84///ymT9OMGTP485//zHe+8x2uuuoqbr31\nVr7xjW/w5je/me9973sMDQ1x9NFH89a3vpVkMtl0u3sDpDSt9hAhuh1jJfZfARcAV1R+T2pTLp/P\ns3jxYjZt2sRRRx3FySefDJjEfv/997NkyRIAMpkMGzZsIJ1O8453vIOenh4AzjjjDFd755xzjv33\n2rVrefbZZ+3/R0ZGSKfTrFixgs997nOcd955nHXWWcydO5fly5fzoQ99iHK5zDvf+U4WL17savfR\nRx/lXe96l02sZ511Fo888ghnnnkmBx10kL39UUcdxaZNm1o6B2eddZa97y9+8Qv7+H/1q1/Zs5FC\nocDmzZs5/PDDW2p7skMikTKM5gnR/Wia2IUQPwVOBGYIIbYAl2ES+v8VQnwY2Ayc3ZZRNWlZtxuW\nxj48PMzpp5/OjTfeyKc+9SmklHzxi1/kn//5n13bX3PNNXXbc1q0hmHwhz/8wX4JWLj00ks57bTT\nuOeeezj22GNZu3YtK1eu5OGHH+buu+/mAx/4AJdccgkf/OAH7X3qWY3xeNz+W1XVhlJM0P6qqtq+\nASklP//5zzn00ENbamtvQyjFhJgsaCUq5lwp5RwpZVRKOVdK+Z9Syl1SyrdIKd9U+e2NmpmUGBgY\n4LrrruOqq66iXC6zevVqvve975HJZAB49dVX2bFjB8cffzz/8z//Q6FQIJPJcPfddwe2ecopp3DD\nDTfY/z/11FMAbNy4kQULFvCFL3yBZcuW8fzzz/Pyyy8za9YsPvKRj/DhD3+YP//5z662Vq5cyV13\n3UUulyObzfLLX/6SE044oQNnwsTq1au5/vrr7RfKX/7yl4711c2QmFZ7iBDdjq4pKdBtWLJkCYsW\nLeKOO+7gAx/4AM899xzHHXccAKlUih/96EcsX76cM888k0WLFnHggQeybNkyBgYGfNu77rrr+MQn\nPsHChQvRNI2VK1dy8803c+2117Ju3TpUVeWII47gbW97G3fccQdXXnkl0WiUVCrFD37wA1dbS5cu\nZc2aNRx99NEAXHTRRSxZsqRl2aVZ/Nu//Ruf+cxnWLhwIVJK5s2bx69//euO9NXNkFIiCaWYEN0P\nMRHOoGXLlknvQhvPPffcpNRsM5kMqVSKXC7HypUrueWWW1i6dOlED2tSYDJdcykl/9iZRQjBQTNe\nX07jEN0DIcSTUspljbYLLfYx4qMf/SjPPvsshUKBCy64ICT1vRSW/RNGxYSYDAiJfYz4yU9+MtFD\nCDEOcNK5lDKsdROiq9E11R1DhOhmOC11IzTaQ3Q5QmIPEaIJeC32ECG6GSGxhwjRBJxcHtJ6iG5H\nSOwhQjQBZ/x6aLCH6HaExO6As2zv2WefXVOjpRlce+21o9pvMpXPffDBBzn99NMnehjjCrfFHjJ7\niO5GSOwOOMv2xmIxbr755pbbGA2x67recvlcXQ8XfBhPuIg95PUQXY6Q2ANwwgkn8PdKhUm/ErnZ\nbJbTTjuNRYsWMX/+fH72s59x3XXXsXXrVk466SROOukkwCygddxxx7F06VLOPvtsuyzBvHnzuPzy\nyzn++OO58847XeVzH3jgAZYsWcKCBQv40Ic+RLFY9N3HiTvvvJP58+ezaNEiVq5cCZhlfE844QSW\nLl3K0qVL+f3vfw+YFveqVat473vfyyGHHMKll17Kj3/8Y44++mgWLFjAxo0bAbOk78UXX8wJJ5zA\nIYcc4pttGlTS929/+xtHH300ixcvZuHChWzYsKF9F2cCEEoxISYTujKOfVemSFEz2tpmPKIwPRVv\nvCHmwhi/+c1vOPXUUwNL5L700kvsu+++dn2Y4eFhBgYGuPrqq1m3bh0zZsxg586dfP3rX2ft2rUk\nk0m+9a1vcfXVV/OVr3wFgEQiwaOPPgrAvffeC5iVE9esWcMDDzzAIYccwgc/+EFuuukmPvOZz9Ts\n48Tll1/Offfdx3777WevsjRr1ix++9vfkkgk2LBhA+eeey5Wxu9f//pXnnvuOaZNm8bBBx/MRRdd\nxOOPP863v/1trr/+eleN94ceeoiNGzdy0kkn2S87C0ElfW+++WY+/elPc95551EqlSb9DCOUYkJM\nJoQWuwNW2d5ly5ZxwAEH8OEPf9hVIjeVStklchcsWMDatWv5whe+wCOPPOJbI+aPf/wjzz77LCtW\nrGDx4sXcfvvtvPzyy/b3zrK+Fl544QUOOuggDjnkEAAuuOACHn744br7AKxYsYI1a9bw3e9+1ybR\ncrnMRz7yERYsWMDZZ5/tKhu8fPly5syZQzwe5w1veAOnnHIKAAsWLHDVnHnve9+Loii86U1v4uCD\nD7YX37Bw//33c8UVV7B48WJOPPFEu6Tvcccdxze/+U2+9a1v8fLLL9dUtZxscIc7TtgwQoRoCl1p\nsTdrWbcblsbuRFDM8iGHHMKTTz7JPffcwxe/+EVOOeUU2xJ37nvyySfz05/+1LcNv4UqGsVIBy1u\ncfPNN7N+/XruvvtuFi9ezFNPPcX111/P7Nmz+etf/4phGCQSCXt7Z3lfRVHs/xVFcS3l582w9P4f\nVNL38MMP55hjjuHuu+9m9erV3Hrrrbz5zW+ue2zdDOd1CXk9RLcjtNgbIKhE7tatW+nt7eX888/n\nX/7lX+zSun19faTTaQCOPfZYHnvsMVu+yOVyvPjii3X7O+yww9i0aZO9zw9/+ENWrVrVcJwbN27k\nmGOO4fLLL2fGjBm88sorDA8PM2fOHBRF4Yc//OGo5JA777wTwzDYuHEjL730Ug2BB5X0femllzj4\n4IP51Kc+xZlnnsnTTz/dct/dhDBBKcRkQlda7N2EoBK59913H5dccgmKohCNRrnpppsAsyjY2972\nNubMmcO6deu47bbbOPfcc20H6Ne//nVbZvFDIpHg+9//PmeffTaaprF8+XIuvvjihuO85JJL2LBh\nA1JK3vKWt7Bo0SI+/vGP8+53v5s777yTk046aVRL2R166KGsWrWK7du3c/PNN7usfggu6fuzn/2M\nH/3oR0SjUfbZZ5+a2cxkQ5igFGIyISzbGyIQa9as4fTTT+c973lPR9qfTNd8OFdmV9Z8OU9Pxhno\njU7wiEK8HtFs2d5QigkRogm4wh1Dmz1ElyOUYkIE4rbbbpvoIXQNjDBBKcQkQlssdiHEZ4UQfxNC\nPCOE+KkQItF4rxAhJg+sGuxCiNBeD9H1GDOxCyH2Az4FLJNSzgdU4H1jbTdEiG6CBBQBgjAqJkT3\no10aewToEUJEgF5ga5vaDRGiKyAlCARChAtthOh+jJnYpZSvAlcBm4FtwLCU8n7vdkKIjwohnhBC\nPDE4ODjWbkOEGFdIJEKY5B46T0N0O9ohxUwF3gEcBOwLJIUQ53u3k1LeIqVcJqVcNnPmzLF22xFY\nZXuPPPJIFi1axNVXX41h1K9Zs2nTpnFZ9/Siiy5ylQTww1133dVwmxCjRIXLhSAMZA/R9WiHFPNW\n4B9SykEpZRn4BfBPbWh33GGVFPjb3/7Gb3/7W+655x7+/d//ve4+40Xst956K0cccUTdbUJi7xwk\nJqkLEfJ6iO5HO4h9M3CsEKJXmIVE3gI814Z26+Kuv7zKiit+x0GX3s2KK37HXX95ta3tz5o1i1tu\nuYUbbrgBKWVgCdxLL72URx55hMWLF3PNNdcEbufEpk2bOOyww7jgggtYuHAh73nPe+wa7kEle088\n8US7MmMqleJLX/oSixYt4thjj2X79u38/ve/51e/+hWXXHIJixcvZuPGjVx33XUcccQRLFy4kPe9\nL/RnjwVSUo2KCZk9hAPD+TJDudJED8MNKeWYf4B/B54HngF+CMTrbX/UUUdJL5599tmaz4Lwyz9v\nkYd9+TfywC/82v457Mu/kb/885am2/BDMpms+WzKlCnytddek9lsVubzeSmllC+++KK0jmHdunXy\ntNNOs7cP2s6Jf/zjHxKQjz76qJRSygsvvFBeeeWVMp/Py7lz58oXXnhBSinlBz7wAXnNNddIKaVc\ntWqV/NOf/iSllBKQv/rVr6SUUl5yySXya1/7mpRSygsuuEDeeeeddj9z5syRhUJBSinlnj17Rnta\nOoZWrvlEY+tQTr66x/zZOpSb6OGE6CJsHcrJV3Znx6Uv4AnZBCe3JSpGSnmZlPIwKeV8KeUHpJTF\ndrQbhCvve4F82V3QKl/WufK+F9rel6yYZ/VK4DrR7Hb7778/K1asAOD888/n0UcfbViy10IsFrOX\npjvqqKNcZXadWLhwIeeddx4/+tGPiETCXLSxwJBVKSaMignhhCG7L2ltUpYU2DqUb+nz0eKll15C\nVVVmzZrFNddcY5fAfeKJJyiV/KdezW7nVw5XNnl3RKNRe39VVV1ldp24++67+cQnPsGTTz7JUUcd\nFbhdiMaQUprhjjR/nUK8PmBayRM9CjcmJbHvO8V/0Yagz0eDwcFBLr74Yj75yU8ihAgsgess0ws0\nXSp38+bN/OEPfwDgpz/9Kccff/yoS/ZacI7FMAxeeeUVTjrpJP7jP/6DoaEhe1m+EK1DVix2RXSf\ndRZiYiFl99UPmpTEfsnqQ+mJqq7PeqIql6w+NGCP5mCtoHTkkUfy1re+lVNOOYXLLrsMgI9//OPc\nfvvtHHvssbz44ot2CdyFCxcSiURYtGgR11xzTeB2Xhx++OHcfvvtLFy4kN27d/Oxj33MVbJ3wYIF\nKIrSVMleC+973/u48sorWbJkCRs2bOD8889nwYIFLFmyhM9+9rNMmTJlTOen2yClZE+2hDFO2ogQ\nmKmnIUI4IGX3yXOTtmzvXX95lSvve4GtQ3n2ndLDJasP5Z1L9mv3UDuCTZs2cfrpp/PMM89M9FAm\nFGMt21so62wdyjO7P0Ey3lkfwuZdOXpipjGRL+kcML23o/2FmDx4eVcW3ZAcPDPV8b6aLds7aT1q\n71yy36Qh8hCdgWWTjIdpYmWeWn+HCGHBstZlpVBcN2BSSjGTHfPmzXvdW+vtgEWwxjjMOs1aMeZP\nt027Q0wsLNWjm+6LSWuxhwhhW+zj8ECZmaei0l8XPcEhJhTOe8H8O7TYQ4QYE2TNHx3sS0oEZlSM\n9X+IEE4rvZss9pDYQ0xaVKfAnX2irH6s6o7mZx3tMsQkgcti7yLfS0jsISYtpOd3x/qxKjsi7Jl2\n9zzCISYS3bpkYtdq7JuHNlPQCm1rLxFJcMCUA9rWnhdf+tKX+MEPfsCePXtaSgTatGkTv//973n/\n+9/v+/2JJ57IVVddxbJl7gin2267jSeeeIIbbrhhTOOezLAepE5b7EaV2auRMV2kp4aYOLgWOe8i\nYu9ai72gFUjFU237aeUl8eCDD7JmzZqWxnvGGWfw+OOPt3iU41f2d6/EODlPreZF1WDvKj01xMTB\nee+FUsxeiGOPPZY5c+bU3eahhx5i8eLFLF68mCVLlpBOp2vK/ubzed73vvexcOFCzjnnHPL5av2b\n73//+xxyyCGsWrWKxx57zP58cHCQd7/73Sxfvpzly5fz2GOPYRgG8+bNY2hoyN7ujW98I9u3b2//\nwU8QrAep045Mq3mlUrbX2XeI1zdklzpPu1aK2Rtx1VVXceONN7JixQoymQyJRIIrrriCq666il//\n+tcAXH311fT29vL000/z9NNPs3TpUgC2bdvGZZddxpNPPsnAwAAnnXQSS5YsAeDTn/40n/3sZzn+\n+OPZvHkzq1ev5rnnnuMd73gHv/zlL7nwwgtZv3498+bNY/bs2RN2/O3GeCUoWSQucEoxHe40xKSA\nW4rpnpsiJHYHjjnmGIrFIplMht27d7N48WIAvvWtb7F69eoxt79ixQo+97nPcd5553HWWWcxd+7c\nmm0efvhhPvWpTwFmHZqFCxcCsH79ek488USsZQXPOeccXnzxRQDWrl3rKg88MjJCOp3mnHPO4fLL\nL+fCCy/kjjvu4JxzzhnzMXQTbOdpp6UYx7J4ItTVQzjQreGOIbE7sH79esDU2G+77TZuu+02+7tX\nXnmFM844A4CLL764peJcFi699FJOO+007rnnHo499ljWrl3ru11QWnLQ54Zh8Ic//IGeHnd1y+OO\nO46///3vDA4Octddd/HlL3+55TF3M8Yr3NGCQIQWewgXpFtk7xqEGnuT2H///Xnqqad46qmnRkXq\nABs3bmTBggV84QtfYNmyZTz//PM1ZX9XrlzJj3/8YwCeeeYZnn76acCcTTz44IPs2rWLcrnMnXfe\nae9zyimnuKJjnnrqKcB8EbzrXe/ic5/7HIcffjjTp08f1bi7FeMe7uh4r4YaewjwhDt20T3RtRZ7\nIpIgU2xf/fBEJNG2tvzw+c9/np/85Cfkcjnmzp3LRRddxFe/+lXXNtdeey3r1q1DVVWOOOII3va2\nt6Eoil32d82aNXzsYx/jwgsvZOHChSxevJijjz4agDlz5vDVr36V4447jjlz5rB06VK71vt1113H\nJz7xCRYuXIimaaxcuZKbb74ZMCWb5cuXu2YfewvscMcOz4GdD6xF7t007Q4xgehSKaYtZXuFEFOA\nW4H5mIf6ISnlH4K2b0fZ3hCTH2O95jvSBTIFjaiqsP+0zpXRzRY1to8U2G9qDwLBlj05ZvbF6UtE\nO9ZniMmBoVyJ3dkSQgj6ExGmp+Id7W+8y/Z+G7hXSvkeIUQMCItVh+g8xjuOHVGtFdPZLkNMEphr\n4Zr3RTdZ7GMmdiFEP7ASWAMgpSwB/gt9hgjRRljP0XhlnpqLWYe1YkJUUS0OJ7pKY2+H8/RgYBD4\nvhDiL0KIW4UQ/uvBNUA3xYGG6Czaca3HLY69WlGgGuwY3qohsCx28+9uoq92EHsEWArcJKVcAmSB\nS70bCSE+KoR4QgjxxODgYE0jiUSCXbt2heT+OoCUkl27dpFIjM2h7cw87eh948o8dfcd4vUNiUQR\nAkURXUXs7dDYtwBbpJTrK///Fz7ELqW8BbgFTOep9/u5c+eyZcsW/Eg/xN6HRCLhm6DVCryV9Tq1\nKpmdeeqQYrpJTw0xcXDO5rrpZT9mYpdSviaEeEUIcaiU8gXgLcCzjfbzIhqNctBBB411OCFeR3DX\nwu5kP+Zvi9SFEOHMMgRg3huKIlCEQO+ie6JdUTH/B/hxJSLmJeDCNrUbIkQg3AWYJGqH0v2dy+KB\nuYpS9zzCISYSEmnXEJJdNI1rC7FLKZ8CGsZWhgjhhW5IippOb2xst2JHJfZK5IMFQXfpqSEmDoYE\nVZgVhLrpnghLCoSYUGQKGq8NF0aVPSolqIqleXfuqTIt9ur/QnSXnhpi4iCltH0vIbGHCFGBpUuO\nhpitiIROQ0qfqo5d9BCHmDiY94b5sh+vYnTNICT2EBMKwyb21ve1HFfOdjoBiayx2LtITg0xgZB2\n5qnoqnd9SOwhJhTGmCx2qin+HQ2Lcf8ruizLMMTEwbCkGLorwTIk9hATCjt7dFQWu0S1l6rrHJzZ\nhVCJiumeZzjEBMI0LoQtCXa60mizCIk9xIRCrzwIrVrAlnVUTRjqtBRTZXZBd027Q0wc7IipLisO\nFxJ7iAnFaDV2i8etqJjOhjtWJR+oRMWEJvvrHpZ17iw10S0O1JDYQ0wo7MUyWnwgrK2rGnuHwx1x\nWuyhFBPCYZ0LbCmmW+6LkNhDTCgsQpdGa/tZRK6Mi8Xujooh1NhD4CnnXPmsW5zqIbGHmFCMWmOv\n/DZjiEVnNXaJK4q922pvhzBR1PRx7c+65ZzO02554YfEHmLC4IwgGK3Gbq1e0/HnSbj/7JYHOISJ\nQlnn1T15ciVt3Pq0LXa6ryZ7SOwhJgxOK7t1jd3xUHW4dos381R0WTJKCJPYoToDHE+EztMQIRxw\nPoMtE7ttsXc+SqUm85QwKqbbkK8Q+3heFdc9SOfzKVpBSOwhJgwuMh/lEyEQlaJcnYPh0djFOETi\nhGgeUkoKZdP73qoTfixwOk+V0GIPEcKEW4ppbV+3xd5p52ltghKE9WK6BUXNsF+y4+nUrjrwRdct\nch4S++sM3ZLyDFViVJXWidn5AHcyxV86HGQWhOL+LsTEwtLXYXyJ1bpnFUe4Y7doMSGxv87wyp4c\nI4XyRA8DqD4YoyJ2j77ZKYvdGdJmoRqzHKIbkC/rxCJKx2duXrgis8ahymgrCIn9dQQpJbohKWvj\nKETWgaWHRhSlZUvLOQ3uqMVu/eEqKdBd0+7XMyx9vSeqjvuShd7ZXDdFS4XE/jqCJX10ixpjT2WV\n1klSOhxXHVrqtLafCroty/D1DEtfT0TVjs7c/GDP5irWutJFi220jdiFEKoQ4i9CiF+3q80Q7YXt\nYOqSm0+X0s7aG22tGCuOvWNSjKMfC92WjPJ6Rr5k6uuJqGpel3G8JhJqq352yT3RTov908BzbWwv\nRJth3XPdZLGPmtgrMo6dedox52m1Hwt2zPIo+zQqC3iHGDsKmqmvq4oY95WtDO8i511U9bMtxC6E\nmAucBtzajvZCdAb1KinuyhTJFscvHdsaj6KMrkKjJYMoHQ53dGa4WrAt9gbm4VCuxNahfM3nI4Uy\nW4cKXUMCkxVOfR3Gv4aPWc7ZmZHcPQ71dlns1wKfBwK9ckKIjwohnhBCPDE4ONimbkO0Auum9yPB\ndEEjVxpfK9Ky2KuLZTS/r7dWjPlZ+x8rZ/SNhWalmJJuUPJxVOuGREpJSe8OJ/ZkhVNfB8tiHr/+\nvVU/xzsqpx7GTOxCiNOBHVLKJ+ttJ6W8RUq5TEq5bObMmWPtNsQoELQMnZQSQ8pxtyB1Q466zoZz\ny7FKI83AXY+9ufRxKf2PyXqBlfXuIIHJCu9iK+PuPKV7l0yMtKGNFcCZQoi3AwmgXwjxIynl+W1o\nO0QbESTFTFS0jLUy0WhKnjqzQa2EIUNKlDaHyDjTxi00W1LArjXvyVy19itrBsTbONjXGbzXZryJ\n1fC5rrphUNbLGNJAIjGkgSENdENHSkl/oh9FdD4YcczELqX8IvBFACHEicC/hKQ+8UgXyiRjETsU\nC6pSjLcCXnV5Oh/r0pB22n67YT0Yo6mzIalq851MGPIbUrP9OV+YquP0WfuVQylmTHDmMlT+GDOx\nG9KozGAN+8ciaN3QMaSBZmhohsarIzmklBQNBc3Q2JXVKOuSvIya45PSNTbd0EnGkijqJCD2EN2H\nsm4wmC4i+6A/EbU/d970hiFrsuX8iHV7ukBEUZjZ137T0pDmNHp0Fnv1ge7kgtZWi04nmT3eBrxs\nlW8wpER1zCSscZa7JTxpksKZ0m/+FjYJN0vOmqG5PvNO+LzkrAilkhSnoBs6EVUlFomREAmK5TJF\nzSAV839WMqVMp05FDdpK7FLKB4EH29lmiNZhr0rkIR4njThlC6OynR8varqkjk981JAVTV9xWbKt\nRcU4p+Bmm20coNVPPSmmwXgbSV/dkgHcjWiGnPfkCgxmCygRFTDYnS0xnC9RJOZqqx45m+UAFCIi\nQkIkWhpjRIGootjSiqB7QolDi30vhBaw3JxTE/arhe5HjIaUGLL9Moz18hGOZcVaeigc23bSeeqf\noNRcf9bU514OAAAgAElEQVT593NWg3luNd0gMg5T84lApyxnUYmkShcMirqBJI6qqMQjMXqiEVKx\n1gh6tPA6T8c7QaoeQmLfC6HrAYTi+Ntv9aKgCI5OVIS0mlRGWcva+VA1a0GPBn4JStb/zWrstcRu\n1sfRDIOyLomo7RlrJ9CInHWpoxu6Lzmbsyq3c9FLzpblrAilZcu5XNYoqhrxiCl9qIoGGDXO6k6h\ntk5/94Q7hsS+F0KraCvem8wpzbj1dvy3r8glnZheOis7NqtZu8dWJduOpvhbxO75uNEqSta5Az8p\nxoy91ooGJd2gh84yezeT89iOyxtHbn3udlZ3DLI23BFqo6AmAiGx74WwNXbP506LNmi9UbdTFdfv\ndkLaFrsYlcUtqaZzV6Wc0Q90OFcmHlXsZBdnP+B+gK3/650X53fezQwJkcoLrdnImMlIzrmSTlmT\nDPR2hmbMWVvr+QXtgjfc0TWucRpDEEJi3wtha+w+EoCFQGJ3OlUdcdjttkKcMcgWubSaeWpLMdZn\nYxjPnlyJZDxSS+xBUgz109f9pS6TnMt6mbJhoEuNkUKRRKxsk7JTd3Zpzs4Dpbss5yBkCjolvYPE\nXiOFWJ93nlqt2ZjistirBpEywcweEvteCNti90orjr+dJKo7tgvS4XVDEmnj/Fa3HwxR+T2KOHZP\nuONo17u0M299iNprFdoJJ1KjqEG+LG3r2UnO+XKJbek8htRJlxV6YmbbujTYNqIxoKmUNINSGTSi\nXUvOY0FJ8z+n7YL0EOh4Vt3060K0w8JoE0Ji3wuhBUkxFavbqf+an1e3MQI/b+8YLRK2iL3VdHAp\nJcJRBxual3JsKaMibZR1nXw5j1AUhgtF01quSBuD6QIjhRJSidqWs0SyY0QjogpyuuprOWs6SAwU\noRBRoqTiZgieZkh6o0X6YlH0iGRYaiRj8YazoaJmEI9MnuiZkmYEZt62C9I5bWN8V7aqOv9rj8sY\nhxlDI4TEvpfBKjAFtZaLNUU0EK7s02Ay99+mHXA6T6H1UDFTW5Z2+nZBK5ItaUQiJdtydurOXlnD\nqTuXNYPtGY2eklK1nq0kFGmgKoqdhGIhGyuhCuiLx3zHV8AgopQqx+YgZFvawZ4BlXVJLBJMBNmi\nzq5MmTlTYkQnSWhkSXOH1nbCmSnxl0LGJTDFcR2r/bu+mlCExN6l0HSD3bkSM1ONrTnXfkZVj6iJ\nY0dWCMsd6aIbElUxyd5N5o5923y3mv1I07knJZpRRi9LsiXdJWsEkfO2oTKxKGR0M7N2e7ZEpqyQ\n1SItyxpFxaA3VqInqpCMuYk6opSJKEZNfY9GJVrr+TCs/SOKk9iD27IqROoGRLs4NNKJoiP5qlNE\nW6NluzT2zsLP0LGf0y5g9pDYuxT5sk6moDHQEyXeQqCzZYkrfs5IWV382SvFVIm9urlLYw94WLyy\nhjNiw4+cLSLfkS5QLEvKwiTiwUwZKaFUqbPhR86qotrk3BMr0htVSMXM7ZOxIr0xhVQ86jvOerAO\nze8IvTquBQWoV5xRBrwUndFA0YoZqzWo8liqfO+t8dPNcJYr7pQ04fV/VPMh2t6Vb9/OPsEZbjnx\n1ykk9i5FUHJLI1j6ejSi+DpPrSgU581f1nViEYFmaBTKRVTFlDfSxRLpYhFdaqAopEvCJmcrzM5s\n1z+czswqNS1di5wjiknkPVGVqGLQFzeTS/LREprE1qIbwsMVZlx5S6fKRjVqxb+bmljHymfSCPbW\nVtdzFb4vS4F5HSKqaFi+t6xZ45t4wmgGhpSUdUk8qlAsG50zYB2RQeAurdxx+ERLdZHvNCT2bkWQ\nTh4Ey1LOl0qU9BKqqlAsa2RKhm0tD2bzlHUdQ2pohk5WV9ENnS17SiTjgmxJktFUUgkVKSXpvEG6\nYCCEgqKqCCXqIuexlB+VUnM/FEJUM6Wa2d8Rxw5jK9lq18rx7Uj60oXSII7dGovqeeNYf1mHHlVE\n3WJgZb3qhGzh9EwoLH09EakQ+zhLMePxAqxfzrnj3TdESOxdBKeskS8XKWglk5hRfGWNkl5md7ZM\nKm5KKRLJnqxOUZPEi4JC2QA1ZpNwplRGoBBRFBShVlKxBb2xCAM9EaTU6I1GSFUEX00rV+pIQ1yN\nEI+073bxOtQaEaUXzsxTwLSgR/lE2UTgs781y/GikcZuSy6KcJVQc0oxYDpQi8VgxnY7IbuAMZpA\nsWweTyKqMJzv3LjN2VT1/04Wg/Pt29190zWExgMhsbcZo9Gc/WSNPVmNbNGgJFV6Y6qvrFEqK5S1\nCEoiSipuXsp8sURMhVhEoAiDVLxaQjQeMYvHqoqgqJkOQc2hyXuzKS3yFARr7KOFGa5YtfhFCxa3\nXXHR8ZnSgGjrt1f57fcd/uqw0uBFYi/UDTi4ucbSi6oKhtTRdP88AfM6jf9CzWNBSTOIqsKOeOoU\n0cnKObYwngGG1cS16mfWn9Y1Lmlm5FZPbPw93iGx+8BFyB6CtqIymiFnC81ozl5Zo1gqIw2dZCxK\nMu5/Y5S0EooQLqvOTCSqpO74hDtGhDsd3irwVdXe3VExFmG2++HUJcRdBnfzcexeOQMqGvuoxxIs\newU6TysvoqAYbSvE2isxeQkhWglzLOkGEbX2Opc0STRi6vSdKMbWCRR1SW9U6bgzsTbzdPxegH4Z\nyd77ZCinUdIlc0Nibx+CyLme5axLHaPyELaDnMc2fqtf/+91Q1KoRB44Q8t0A+LRAKKUEhTFZW3a\nWrAiahyQFml14mHxLgTcSgElv3MihBh11Egji91Xi3F87/et9VL0SkxeKSZWsdJLmqTX4zc2yw8Y\npOIqZV3WjcJxQjckrw2XmJGKEo+Ob9x7WTcwDEksonRUGvHTuKE24qtTcDrB7b4dUoxuSFMKpRpO\nPJ6YlMRe1svkyrmG5AzUPnWSQHKOiui4rEfYDIIqA1ooVJxSvTGVXEmvLAxtbh9RPC8k64ajSjaW\ntel8QLwOSAM34bbz2KwQSwutVOazRlIzDR+1xl4dl99Y/aItGtUFsbIivRJTbUVCQTSi+C66UdbN\n8xSLKOiG4cpRqIdC2ZQASrox7sRuzR5jEWHfd50sIlfzzh2DE3008F57UXkGLVIHM1ItJPYmUNAK\nvDLyCj2RHqB7yXksaPQw5Io6EVWQSpjEXtIMW6O1YtLBbVFaU1frYdCluy66EG5rx5JuwK0TjxWG\nz0PZUtagn745hllFvYVGgqSYRgWnrJeid+ZU4/TFtNrz5VrSdpJkUQOrFlgjWM7LiYiisXwCsYhl\nOHWmRrl35mNhLL6W0fQfVKc/XzLsl7qmS+LjzLRjZkAhxP5CiHVCiOeEEH8TQny6HQNrhKgSpSfa\nQ0+0h3gkTlSNoirqXkHqUD/c0ZJheqKKnYpe0iRWBVhVcZTCdcoAgLk0mEWiVYtdEaJSbsA1CHN2\no7R3ehv0UEJzeqxfTZhGUSr14JZK3K1Y58yLRskwUpoPl3fGI6l9UcQiAsOQtiPbQkk3KolMCqoQ\nTWvsljQ3EVE0ui5RHVOuVpzirSCoLpCg+Rf8WO5p78zL2b+U5qzJDHqgYZ5CJ9AOFtSA/0dKeThw\nLPAJIcQRbWj3dY1qLfTam8KWYeKq+eBHFIoVDzyAqjhqUzt3r1jsTlKqkmxt4pIurc/bO522nJWj\nrbMRFJEw2gfVlYXr853fJLpRaJvTPwHu2va1Frv5GHrlmJJWrSFjz7IaXAizfo4l49Xd1IU92TL5\nkt78DgFwynfQehhr0/047lsnRAtvkq1DJdKFJqdBHkhqr6M1nnzZzD3oiSlEFNEws7gTGDOxSym3\nSSn/XPk7DTwH7DfWdltFUTPYsqfo0rbageGcxkh+dBd/LKjn0MsVdVRF2NX+4qqoWOwWsfsvXmEt\nDOCUPXRD2tEbpsbulA1kRXvvvMXeihRjbeJ8rMwxjm48hhGcXGJKMbUPsFuK8WnTdp66j8vvRRF1\nzLqq/ZqOU0vSsDTaRla4df+bwTjNnRBDStIFnWydePpmYXjOl1feaxsCpRCaWnpdN6TrJdhy9zIg\nvLIysxLCjOOPNEhA6xTaqlsIIeYBS4D17Wy3ETRDMpgum9PZNr8dsyWdoZxGrg3WTCswApjdkKYM\n0xurXrp4VLEdNkqFuL1E5Yz9doahOYnLO421vnM6W9t5bKOus+HzUDci2kbjsQnY2U1A5AU0J8U4\n93OWLfBamVbdmJLj3i1VHKdWPZlm66A4Ne5madp6ZrxS0GhgxpZX/++cxe5/bRSaMw4sI6jJBaxq\nEOR7sT5LRMzos4gq0EbbyRjQNmIXQqSAnwOfkVKO+Hz/USHEE0KIJwYHB9vVLVKapF7VMNt7F1k3\nwO6M1vQyZmNFUIVFqCxeIHElPVjT9aJmYFV1tR2mnt+WBQ7mzekkGuGw2N3RMuYG7Xpn+lnsQRaz\n7/7WPo7PvJJHs/CWD/Y6j4PQSIqpzo7cYw6awkdVxVU4q1Ay/7aiWuzlChscYLFsEI+YL/dmLXbL\nam2HUeSc/UDVmdhu+N0DlQ6b6q+e1NlU/943t2c81vMZUc2ZZDtemq2gLcQuhIhikvqPpZS/8NtG\nSnmLlHKZlHLZzJkz29EtALsyGmXNYFoyWumnbU3bFm1fQgUBOzPaOMXIVv/203zBfU9FVcWWIiyC\nqpUAsPdzWsfOdRudcoZzIYGxWMN+sF6WXgJwjrMe7HG4NPZai7sZWJEjqs+LxfrTV4qx9g8IkTT3\n85di/Cy9WETY8oAhJemiTk9MtUNX1SZCBy1pIR5VUFuwlKuavBx1LoAF5+wHrBDa9j8zQeGOzfZn\nHWejl1+upLvyROz+CYiWqvxOVGbUVv388dbZ2xEVI4D/BJ6TUl499iE1j6JmkCvpTOmN2NmZ7byH\nrIsfiyjMSEUpawZ7cp3X2103pueAgpxGltVux8v6RGOYHzusSOmOG3dLNE6L3d33WCF9jkFx9N1w\n/8pvP4u9VRKpsdh9Oqonxfi36djPc1xGgDZraellXZIp6BiGpL+nOitrxnlq6evxStZns9aok3TG\nQkBBEltn49g9GjvNcYAtxTTYdk9WI52vlWGDrmM8qtDreCFHmizN3G60w2JfAXwAeLMQ4qnKz9vb\n0G5DWPG6FqmPJeTND87wwURUIRFVxiV0yVXm1fNd1Vh131Zxj5PNK8U4ScqZTepMYHIm3ThfIPbn\nbXpC9YpV51vytCmLvXos9v4tSDl+bSlezYT6xFgvI9cpNdVkX/rP4G0tvVg2SBd0ElHFtRRe1Xka\nfCy2vl6p09KsX6SsG9Xa8GO4xs5ZnoV2O96rfVUMD+8Xorl7wFlSo1HNn6BZmd9Mrr8nwoy+6poA\nkUogw3jJuHa/Y21ASvkoAQ7iTqNYdhcbanXdzEZwhg/C2NLWvfBOWZ1wEkNNlAa1MgY4LXZc31d9\nsO6GREWDddZscVq9LnJSrDE3d2yNYBigeEwKiyhb6cJdi9tEq0PUpfsaO/evSjFw7zPbuOnBjbw2\nXGCfgQQfO/ENzN9vmu/95vRPeCWmoOuuKqajLV3QMaRkeqp2wRCztnvwEVr6ujPMUpfVJDM/mKtX\nyUrZAn1Mhouz7pA9ZlF1vHdi3VM/R3Qzfjbncxx0jiwp1u++D8o49kNEVdAMaH0JmNFj0mbzSCnN\n8rSOlGkhaKvJ7gwfhIqHvwG7GbI20cSLTEFn61CpTqhcZUqr+Gns2GNxIh5V6ImpJCrno5ZQqHxe\n3d+y2C1r1fkyqEdOY0UQuTVbCMzexEUgoxuj15HrlsHMX2uf3c4373mebcMFJLBtuMA373meR17c\n4dufn8XutPyC+M2s9CiJRhT7OjqhiuBsUqe+Hng8PtCMatkCM4Jj9BfZ6sovfT7okSiWDXamyy33\nVb2f/aWYRsftfI6DnmnrXAe9vJt9T5mx7JPMYp8olHXzQXFOV9svxVTjuO32G3QwlNPIlwz2mxoP\n3KZUKZSkGdKeAjvhJAbvTRV0QytCMNMxBaxasLbJbm9n7V8Nd3R/Zz4Y1c9a0b+bQZCF1OxiGfas\nxfnhGDX2iI/UYX333UdeolB266yFss6P12/m5CP2CWxTCLfEFOTwsxCPCPIl6E/4VwOs9+Jz6uvg\n8FkYQJ3igpaFHlFFw0U/GsHvxWXdU/c+s41bHn7JNeM5df4cskW9Uuso0lI9lSApxDZOqC8jaLL6\nfAUdsjWb837vV+vIhmEgsjlEJoPIZFDSWZJDWdLDeXpJw3kfaXxwbcCkJfZiJZnDZbHT5gxJw33x\nmql7oRvmC6GoGa6Xjncb67ff4sROh57uKdISdEN7USvFVD6v/FYFrlrszn281qVTe28HDEOi+Fik\nzV4/P4IcrYPXssq80hBUz9n2kYLvvoPpou+LxH754paYvLMDL3orviK/Ms33PrONG373EoPpIvsM\nxGxitJArGa6kNTs8spHFXiH2qCpQ6yz64SdFOft3HrfXefrIhkFufujvFDXz5WjNeAAW7z/dHqfa\ngqIrA5jbdd/Xac4wzHLIxbJPJFCFnOWeEdSdGZRsjhg5m6iNdIaetE5fdpj+zB6TyNMZlEwGkc/X\nWCdCjaNFkyRlLiT2RihpkrgqXJUMW0knbgaGlB5ir34eSK6V/k29M4jYzW00XfoKb/aU1udwGt2w\n1bG6pQlviKDpL6hmKZr7mL+9pQagsb5rwbKAInWsL/P81X7e6vVzNtFILgoiJuta+u1v/T27P8Fr\nPuQ+sy9eNwnIOztqJMVEFEF/T+0jee8z2/jmPc9XrHLhIsZT588xZcmyO2mt2RddWTfvcStJygp5\ndN731f5ridlJ7s66Q85x/GT9ZpvULRTKOt9Zt5Ebz5tq7ttgZuGFBN97yL7vdQORqZBxJovIZCtW\ndBYlnSE2XKYvlyabLdKbGWIgM2R+5yDntBqnFE0CkCzstnXrYjROdNocYnEQvXGMaVORB8xFplIY\nqSQylTT/7ksik0l++2qe/++pXWwu5JjzH+v4wqlH8s4lnU3On7TEXtQk/XE3cXZCiomozoelarkG\n8Zb1oOfLBv09we1CcASCWU/EJFOvIyjohvaD83x4LXan7KE4yB6qYZBOfb1ZmWQwXaakGUxPRen1\nWWCg3jS2pghZAPzkKEvR8nNu1yOmo+fNcEsmjvNtHe/Fqw7mintfcMkxiajKB487MMCx5ibwqu5b\n/b8V3PTgxkrf1T0LZZ2bHtzIqfPnUHDUJrFQvVfrX7SyXpUDrZexpruJvdp/Fc7+Lfj5fxQh2Jkp\n+va9faRg7+O7Qpdh2ITsJef4njxKNk+yMISSztifxzJFtJJk6vB2YtK8m3ZnS7y6J09JN4ipCvtM\n7SU2ez964wqF/qnIRARj+jTkgfu7yHkk0Ucm1suTgwXufmGQTUVBauYU1px0KAvnTiPeF6XcYBGN\ne5/Zxjee2k5Ri2BEomwdzvPFX/wvQEfJfVISe0kzNWpvrelOSDFxh0XtjpX2fzyt/ouVh81r2Usp\n3Ra7bxuVAlLUkmlQ/KwfnIsOeOULl4zhCZH0Ji6Z24uG5zZb1CmWTUlgZ7rMlF5ZY4Hq9sPv7zxt\nJurIu0iHta8ihC9B1COmoy6YgVXozGzb0U+F5FfPn4OiiBqLf9mBMyj5OMVqZjuiMgsi+Njr4bVh\nfynI+jxfMsMcnQ5Xl8ZeB5ou6Y17kmkMidND1Kh/C866QxaENDggKins3kOqlCdZKpIqFUiW8syK\nQs/Pd0EuT29mN/3ZEduiFpblHIBsrI9INEYsbiCTSYxUCmPGNLQD+yn1TCHbH6Xcn+TxnSX+86kd\n7FFjZGIJMrEe9J5ePvSWw3nnkn3J5jSMqILiE4WUz5ZZ97+vcfOmlyiq/dALu7M637r3BS5e9Ube\nc1RjYr7pwY2V2YqC9YTlyzpX3vdCSOxe2I6iSC2xt4vXrZK2EV8ppv5+Vj10q3SnE04uD7LYLYem\n80VStaZbCBtzWNlVa9GywGuPy2lpOZ2q1nf1rD9DSoZyGrGIwqz+KLszGkM5DUPClN7qbWZFIPhp\n2oL6zs+qnFJkZirBJ99ysMtiVBV/IqtHTFZUkNPJiedvIUzJwasp786UfWcx3tmOEIJ1z2/nR3/c\nxPaRErP7Y3z8pFqNOgj7DCTY5nMM+wwkAHN2aCYluV/EjXxCWiXLNWrX8acSc+3eZ05fjJGdwxVi\nLtjkvH9Up+dnI4isaU0XR4qIXJGB9KBNzlqhyNVGkuy2HSTK1WNQhKBvn5moI5swkr2QUKA3hjFj\nGvq8A0yy7jMlDdOK7jX/7uvDSPaymzgyGiXR7152qlA2yI+UyPfHMKIK/3bDo2zbp6/m2H+yfjPv\nPmo/FCXYmDCkv4xU1HR+sn4z713emJiD7r2tQ8EvrXZgUhJ7UdPtuF8nvOtLjgXO5CQLzYTUGRJ6\nowq5kkGh5EPsNrGJBlKMu+SrnfJO81KM4pJi3PKA28FVJSAhqglKrvUcqb/YxnBOQzckM/qiKEIw\noy/KYBrSBc1N7Lb/wF+KqaeRV+UUwY5MsUbnDfID1CNG69w6V5my0MjZGeQSMKR71aWHX9zB9b+z\nnIeC10b8NeogfOzEN7iOHUwp6GMnvsEu1+zU153js873fX/dwu2/fZbs4BAHxiUfnD+dpTMSxPcU\n6Cul6c2lEZksA2mNeDbNQGaPTc7/vXOYl3flXOdWEYIDp/fSuymGjMaQfUmUvhmoqT4XOZdTSWKx\nAV4dLnL3izvZXFZITJ/C+958GEccegDEopQ0AxlXIdl8pLcxXPKV5KuSmokgYt2ZKaIq5rUPivDU\nDRkgI5nyUjMGVvXec8/y950SoNO2CZOS2Atlg3jUbyrfPinGJmAXAZq/61lBlkwQjyq+q+JY7cYj\nCvmS7muBW/VE/OqfGLJ+wokTTimnnhSjekjeKing7EZRBNKnZgaY0li6oJOKq65ZlBW+55Sk9HoW\nex1ib0bnVYT/y9JNjCYsYjRnZdXYb29t9nrPbpBFLKV7vx/84eWaeiN+GnUQrG1ufmADI7szHBTT\nuOiImawsbmf4wY3EhwtMLY8QyWZcUkY2D9FMGn3LFg7fuptvOon5l4Lc7Gn0DExhoLAHNRpF9iWJ\n98+knEy5yDnRlyS7u8wv/z7MlpJCYvoA71p1GAcunceuvhTETKt5z0gJQ0J6oGpFSykp7i5yRG+E\nYz2y3JbdRXpVgaY3lvm8zu9zj57HWw6fXbOd9xn1f6kLZqTiKMJMbtQC7mtdmg7yHekiXj1gRire\ncAlHcN57hr1/T1TlktWHNt55DJh0xF6qWChxH3ZrpxRjEUSrFrspYQh6ooJ8Sa+ka1dZrErsJun5\nxbJXrUirTcfbXkqEHyv6QdSmc1edp5aV7rHMRdV52mzNj3RBRxGCgV737WSdO92QKKrbmddqglIz\nOq8qoOQzSIsY/aJituwpVl8ywmux+693Wh2vYzvH8XhlrMF0AevMK4akt5yjr1gguatA9MlUjXPQ\nJGZHFEcmw/szWc4ulHg1PsCsUprk782knnS8n6lS0ldKI2NxZEWyMFIpmD6N0v77clehl61TI2Rj\nCbKxOJlYD9lYgtj0qfz7eccyPHeqTc7ZbJlc0SA1zZ2HcXjlxwnvmfYLKgha97SsG3YeSlELjiUH\nf+f3jes2ElHgPcvmevqrXgPwf6nHIyrvP+aASjRQcP0dw5B88LgDuf53GylqWs3+zVjszntv63Ce\nfQd6wqgYP0gkvTGVvOE/9WwXs3uLQ9ntE2yxO6MhLEdWvmQQ7XESu/m9ZdlquqyJZZcV8vZz6LWS\nyuw8HVWL3dLY3f8796lmpLqjLIIy9MzoIVET6RKxib0ar2/HjfscgxJAlNBYZwbq6qV+GrnVl7Me\nvfP+CfSR6zoik0XdNYK6J4f6jyyRbNYm5/hwgXgmR19uDyKT4ZtPbaNcKDM9N4KmRpiT3gmYKyf1\nvzDgHo+HnI1ZM9APOhAjlaKcTJKPDzDSn0CfkqScSrFLJuifnmLXjAGIuqWMXLqMphtcf+Mjvo+F\nQCEyvUrqgB3yqBmybsiqH7z3jAW/mY21qEhUFZVKlMEPrt9srajp3P77TT7E7ja+/F7qHzh2Hse9\nwawwa9XU8Qt0MKTk5CNmU9Yldzz+MjvS1f2PfUPzFWqtey9TynDQlIOIqp0vLjDpiD0eUZnVH2dr\nujWLr1VYBOwXxx7UhbMIUkQ144ILnrBH3TBvIGtdSD/pwCJvr15o9z0qjd0/I9DPwpKVWHT39mai\nyc9+8jLbR7yx4P5E7bTYwbS8bvzdRnakS75JNn4+BQv1dGbnGK3+GmYxahoynUG8NkJULxAtZont\nzBIr5OkpjCAyGQojJdRcgf7MzmqMcyZrR2vE1Bh6NEWqOERUVqfzseQMSCZQEwIjleLAQ+Zy/2sa\nzwJDiT4EOuWeXt594qEcsHheJaojiexL1ZCzE0ZF1ij0Roj1RMgUdPRsmcSUWLXYjfP8CyjJ4Jfi\nzL54zWzRfhnrrRN70H3gJ7GVNHOx56hqOnnrFckKmq2ZEomnr8pvZ9iq96U+OFKy/UV2cTUDFIeB\nZWVlq4pg1SGzePuCfZha8QHsTJfHvahXq5h0xF4P9Sy+VmERsBOWczHIKvTWSo9HFbJFt6VhkY5V\n9c0v5NHSaH2TZmgu8xTc4Z9eh17Viepuq1pDxk3sv3t+Ozc/9BJFzZQAnLHgi/af7nu+La7RDdkw\nyQaC46id25iWV4mZfXH+z8oDOHVuL+LVrSiZLIk9GaJ7ssS0LLFspipxOOKclWwWkc4iCnk0BOnE\nVKaWs/TpRUbiA0QNnd5yBhmLow7MIJLqgx7hspxlhYRziSSZSJKhWSmi/SmbnHfnTFKMV0o8zMhp\nDDy9jd88/rLrpXbU/Dm0UgTaqo1vXdN8JbQ06kPqYMlnMlCOOPfoA2rI22rLrDvTwuAIrp+iCD8p\nRhJVFazFSOqVqAl6Mc3qqy3b0cj4AndJC9sYkJKI4/kwHMET3hdTpwqatRN7FbFbqGfUZiolUb0R\nNQO0Y58AACAASURBVF4YAVafadEG9OtxUMZUhYx06+y6Ie2/I4rAz29TzYa02pX2bylbSHARAmnd\noZ6T0thid79A/vPRf/hmD9704EZuOn9aoGZuRf9Up9OKa/9bHniRt83tRWSz9KYzRLZnMUoZEiUz\nTVuk0yiVlO33ZjKck8myvQhKJsec+9wrccWUCIVYP4nSCAlDM2WNSsickUphzJ6JnjrIJudSb5Jc\npI/MtF6MqSkyJCDZS2TOFIhGGapYdsmBGH4olXS0dJlyf8xVIkHKouslqgg44U0zWX3kPpR0g32n\nBNcRagRLEvPLNvXCkhlWH2nWs3HKEf+88mDm7zet5h6vhjy2lgpaL/HMb93TklaN5FGU+qs9Bb2Y\nPnT8Qb59meMJHqtuSKINVqay8iEUQc2Lx1rkvZuxVxG7S5P2OfGaIdmdLdOXUO1pVRC0gOgTP+vD\ngjc8zlqkuOzQ0XUDEpWuIz4hj74FpDy/mzUWvFKMX4RPbTEx/4igHSMF8CkG+tpwAaOsoeSKKMMF\nl4WsZHMkdxeIZtO88/4nSRYLaEqUvmKeOZmdpEoFElqJafdPtdvLVqzmZDljjtuHnMu90yDZS24g\n7kjfTlHs6SUteojO7KNnal9dWQPMJLLSSAm9L4oWUzFGSuY5r+znPWe151fY2znh53gGkyzGaukp\nikCX5pqohqxN0nOi2m+tHFHUDLYPl2oUHCHMejOZok6vJ8qpHgyPUeMas3CfI60ydmtxkUYSmlcn\nn92f4JzlB/pGxVhjaJRzYfVjhd16lRVnWLKieKKlAnwJ3YS9itgbRa38+q9b+e7Dm9iZKbDPQNy3\nkJEF3SezFfytDwveVPJYZVZQ1iTEqhmd9k2lCoolw9NG9Vi8zlPD8+JoBOHY12vpW8+PPXHRNNMh\nODiCGMoSyeaIa1kSBZOoL37+fxkqSmal9zC1kCFVKtBbLjDdKDHw55/RrxWZquVqxpCK9aHH4hyz\ncwc7RYwtA7PY3ZtidzJJJpYgMtDHYe9cYpNzSe1lJNFLdE4/MpX0JefMcImIIsj3ub8zDIm+p4ie\njPpXV/OgWh5ZVM+X6/v6xB7kTPdmxjoXKhkzsVdIy1oT1a+8r7dfv0xpJ3F5MT0VZftIicF0mX36\nYw1nt9A42skpX1qhn1HPql/1CoE5X0yaIdm6p1gnvyB43VNvRVPrxVZTRdUhxSjUWuzxsV3GjmOv\nInYL0keMufeZbfy/92+o3FT+Gq+9v5SBUowqguuZeKvbCeFeff43/7uN76z7BzszefYZSLDmnw7i\nqAOnu2QPp07vnVb6rXfqQoWcLW05ujNDOZ0joWeI78kTyeZI5YcqFnWW4YJgSnqIacO7EMWKhhnp\nRYmYkSb95QxJvQTA+8sKD2l9CK2MpirsSA1QTOzDiqXzyM87kGhfnMyUJEZfqmpFp5KMEKOAyuZX\nd1em0xKTPiWJqMq/vv0wCo7zr+Q1SjmN8kA82HkXEIaoBhBtEGquF+7ZWCM/dZAhYYW8eqEbEBvj\nE6cIs7Z3QTOIRZS6TmJXWQHPe87wScCzoCpmCejtw2UG0yVmD8QaGhN+lR2d43AaQ1Zmq2X4NJPR\n7UKd2QHURjc54V1jwS5FETBzVoX58is7NFPnGgbdir2K2P3CAy2YNRuaSxKx3s5BemGw87S6jYWI\nqqDpBvc+s40rfvOCXW5423CBq3/7IheveiNnL9vPXgHJKGuIdJrISJFoIUtke5ZIxXJWh9P0jJTp\nyw+TzJh1NUS2Eq2RzlbJ2TqOSC9E4iQLe8jE+tihKazfNcQuEYNUkjccsh9zFhxCoVKFTqaSZBMp\nMrEUsreH+IwUpal9yFQSVY0w+PgW7vjTy+xwRMUccOQ+bNldJNEboehTmVDNaeh5za5CeMPv/lGZ\nMfmXfrVmSSXNIBJQYMkb4eO8NkFJSn6ozrDsN3HVJwENcwb8LHb/muTV74QY2xReFVA0TA3crxKk\nE/VK99oEF8BPUVWpZA+XGMpp9mLxQfDKkE54cyA03QyPtUNvA3TuIDR6cdeTS/1eaIriL8VYVT+d\n45cei79bsVcRez0pxgyZ8teIwSxgFVFNfTEoO9IM13uJHQE1se2HWtMQObNkaM+uDCN7svzprvUc\nP5wlYkhmZIcYKOaIl0vE/pCgd26EvvQQIp2lUNYoxvqZVhohbmgMJ6aR1HIktQJb8wZ/N1JsKOSJ\n9sR548GzOeDAWehvPNiMyEiaURlGXwqZTJKOJRmOJNi93xR+8dwern94kyvRIhFV+ddTDnMdg5bX\n0CoLdsv+GNJyMknJCW+ayRmL5rgIxVvT3QtrqqsZkpOPmMOR+05jWjJKKmAxiZhqPkimc813k7oO\nZCWgXkxQO1AlN68W7A359CIoz8Bsy+k8rf49Vj5wlk3oqSPDmP26x+SEblRq5NQ5QGvd1VK9WhKO\n9iBIY3fLl2VDEnUwo6VzN2uxW5sFkqtDipFSMpguM9ATIR5VHPere3x+C9pYz78zh8My+ro8KKY9\nxC6EOBX4NuaE71Yp5RXtaLflcVR++90fs/oSbE+X8E6wrQSX3dkyEUWwz0DMvIi6TmQ4h1rIIzIZ\n1j+9mXUPPs8RRY0leplZmT3sXFtkx/5J3tgjzEJIuTKyIJk5ss3uIaZEMWJ9nLPhZUqRKOl4kqm5\nYbKxBOlYL1vUWZRm7Uv5DQdgpJJke1LkE/1kp/VQ6kuSphcxNclvB/Nc/uBmcoaCJQZZUkaQn8DI\naxg5DWMgzg/+9ExTMxYXCfnpxDWyQ315SHVYY/YDWYePTPlKcY21WDbYmSkzu6L3SghkSDWgXowf\naop14Q0tbV2K8T0fPtb7aKHYLyFhz/Iajc/XYpcBKwB5EFEFuYDFN5ywLXafNq1sZitMUNclcUfJ\nbet+aHY9Yb+ZsRPmdTQ3KumyEmKrMSsa8008VH2CGJzhzs4Zl9VuqxU6xxtjJnYhhArcCJwMbAH+\nJIT4lZTy2bG23fpgzF+2daBplZVNsnx8X8Fdj71CIjtCxJBMzw0xQytw+rw+ei5/jGROoOSyxEZ2\nEsvkKMsoMwt7iFToaM6WYT6mG4zEesnGepmWHyYTS/DycIo3Lj8YY/YsSqkplHr6yA9Ekf19GKkk\n5Z4kaaWXb933ApuKCrlYD7qD2Galenj7RUejVKa6uaJOMVNGmxIHVWDsLqDFVb5973pynuer2Xoj\nhjRX+/F75XmTP/wkBPt/H9Ksp61CbZIS+BcAcyIWEWSL1To6w3mzwFi+rNOnRqCexd6KFON1Znos\ndhpY7FZ/bimm0lTAeRwrH1iEYi1aXX/byph8eFk3ZFO1TprNRK06T4PbMiQg3RUlzX1EjVxTD42M\nCefMq1yZbRTKBpoufYv7+ZWicMe6V42aRvd7t6AdFvvRwN+llC8BCCHuAN4BdIbYh4ZQf/8wqcHN\n9BR0R0Zghp5MnnwBBoZ3khoZQpTMzDQDWJGYxj+NjDC0c4itiQFmlbLMm9HDjB3TyWX7oW8GxowZ\n7Dxwf3pScQqJPgrTEmaoXTLJ53/yv6RjCTKxXhc5C2D9l94KQDpbplgyyHvWOzV2F3j72/q4du0G\ndM3AaXGff9yB6LqTGNw3reXhNwm49m4KysoDdzieWcyofkq+uU/1b79QSD9Hod+2FqrE7kiKaiAz\nxyMKmYJOWTctJKtMc6Es6UuYzvHA/nwe0iB4M1y9TjdvTXo/eJNXvGvLmtv4yzKjgXXuEnXi1539\nBhUqc8Zy14Ndo71BJmo9jd05s7Fe8N62hKgfy+6E/fIM+F5QlUzKumFfo0xRN5dl9PTtV4rCuXyj\n7UsxpKPsdHczezuIfT/gFcf/W4BjvBsJIT4KfBTggAMOGHVnD6z9M1M+/2WKepGYGmH2PlOYMWc6\nRjIJqRT6zJmUDtqXQl+iUr85RSbRSyaSYtbMJAP9SRQtxsCMftRUnD2Y+nomU2ZKb4ShnEZGMcm0\n6CDo9PpcJftNVH5MsnESY5D2G1UVTnjTTDRd8sM/vsxgOm87D486YLqrHK70WAQW2ZjZd7Up1F5i\ndsJ+xiScd8wBfOfBjTUauzMl39zHSULe9mofvoZSTOVzaxEGs936D4UVO13SJPmSjqIIEhGFQtlo\nmKTV7BJ+1thrFxORruNq9Px6Lc16zlPv56NBLKKQiCr0NFi5J2h8FrzLPgbBmQ3st4yjBV0GZ0Q7\nnczWbMobQlmvfG4Q6oU7WtMUK8NVVQSZok7cp6aRXxy9KyxZqb6YXk8Wu98h1lwiKeUtwC0Ay5Yt\na/ESmrjrL6/ylSczJE79OCOxPPnYFKLxmK0za4Yku6dIPBkFh3NuJFPGKBmoU2MgBGJ3Ed0RnWCF\nX/UlVLIlg7JmEPUkZrjLb5rwEmNQ3HNMFeTKBse/aSZvPWI20x2rtezOll0LCNfovpUH82MnvoFv\n3P2CS3v2I2Y/GFKy8pBZRFXBD/+4qe6CxH66ugW/UM9GN7oQwl54RPE8KEGIqGZSSKaoU9IMpvSa\nK9jnSrrtyKs7DZfN1YvxZhDa70HpXZAwGN6Yab/z4ZLbx0gIEUUwqz/Aq+wDv6xO3Wg+siNScWY3\nykSVMngmZvUjqT5rXou9tRdy5Y/Ae84hxeiSRFShN2aWyc4bsib23xtH7411d76YaiKpuhTtIPYt\nwP6O/+cCW9vQbg2uvO8FRqTKUP90NHajoKI7dGbrVHtvkILmXmFGVXCl8jvDrwZ6VHamjRr90SLA\n76zbWFkFJ1G7Ck7A1D2iCoxirdMGzBvcqWH6LUknK/3nSwa3PrLJjoOvl2AF7gdKSnjz4bNrquEF\n7dPs0nWNpBjAJnZofgobVwX5srnkWyqhuuqj1EMzyS7VsbszCIXnfDk/C4I3RtvvfFghc0Hx7Z2E\n38vYG8vdCGbIbn3SNQ0S//acUow31NG5jdZkYa1Ger6ojEc3ZKWEh6AnphJRNd86RM5CYKi1IZHO\nqB1vJFW3oh3E/ifgTUKIg4BXgfcB729DuzUIWk7K0pn9LrQhJZouSTlSxZxEA7gcQ70xlURU983o\nO3X+HE48dDY7RkrM6o/VbGMEjCHmsP69lor1XbkSt+21pEwpxhzrSYfNZtm86ew/LVh+cUI4mEoi\nmwq1c84UvFBEbaxxM04zVZglGoRs/oGIRcyFSlIJFUWYhaKiFbIPGh/UPqT1EJT6LyV1NWMnvFmJ\nQeejXsZyJ2FGofhnNzdbvTGqCN+1Xd1tyjqztuo23lBHC6qAYounJ1BjF4CsFtizHLWpuMpQTqsN\nCqj8b6316i2p4bbY3TPqbsWYCx5IKTXgk8B9wHPA/5VS/m2s7fohaDkpS2f2K+rvF1ccqUPsALP6\nY4HJH86L7EWQ1eKMAPDW5rBC1iyJwduGcFhc3gqNjWBLCzRvLVYtdv/vAp2nddq0XqRGCxl7vXGF\nnqhCn0NSi0cVOwMw6DzYemkDEvUbj3PFqqbFGOGtI2K15dnMM77xgipqNfZ6q1j5IaIKNF3WfTEZ\nMphMXM5TXdolq13bKM2/+BqRq+UrsV5GlqyajJtGgrcapjeRyy4AVvncGV3UjrIQ44G2VLKRUt4j\npTxESvkGKeU32tGmHy5ZfSg9nhogXp1Z1EyNazVZi2ikNKdqhiGbqocBDnnD5x6UAVaLqogaR0y1\nPXfZAW8birMv2aJGa5FcneSRgF18b17/RI7alem9UCs6r18p5CBEVYWZ/THX+XLOkAI1dqfFXgfD\neQ0pcVVHdDqbW5FivM5Tv/NRT+LqJEwfTa3GDs1LMZZhUi+MtF6avXXIVvEvv5mCZTQ0o7M3WtnK\nar6kmfdbxPHs7Ts19v+3d64xklzVHf+fW9VdPT0znt31erFZbLAJSgQikckKJQQFCASME7Dz+ECU\nSLwiy3kpUZQXWSmK8oVE+UAUESlyEFKQEJCQF4mChAmOIgUtyEHYhhDANoQQEzDY693ZmZ2ux8mH\ne2/Vreqq6qquR9fU3p806p56nr59+9S5555z7lJyXFyKIkq/ZssOaIt96G4YoCXF3hd3334W7/zJ\nF+PZe1sgADftzZYSdKpMZpkheEWz9EWUJX1EJdEa+seR1/knjsBCWaJLlo/xeVZlQi7LKl/reAC0\nPzjvN5o3WmGsHpbq9vbDarHTRZQVvIrvVTKi0iyCCPtXQ+zMnJSbTMPgyq4Yc6JOnpvfHmUuri7R\npXtNY0dHKFV9yOjfhl/iZ4+ikodtrGiVBZ3TCZySB7IZ8ir/r9aOiyCKC40lsuR/NxNX4ECtQWzW\niUmOSR48x8FiP3YlBe6+/Sxe/cLr8MTlJ7Dr7S7tJ2RSvHMm7XRHDSMj/Kqi9bLKYi/6zqcqNTvv\nPlOXcLDgeBRh1ichJD/KugX+tVWTF4JXeh7lJ7/krTZTJbrCHPk2if+V2ZbyIVgWhbMqSeniQQAh\nCHsZd1s85OZEWa+SNpsuH0XFbiwp34oLtoy+n5lwE0bVrXXAsNjLFHtpuKNyN4bFRpTpQnMzrX7o\nR/j2ZR/P2pvCc0WlRcYB+SDaqWAMAMCJLQdPXvaxfxSm6sSY8oW8nPswVI6VxV4F/WTV5Fvs8jVQ\n7gEg34rIv35xllyZH/u6LQfP2ssPBJ7GcdvR0oSe6ddmrFaiaWHlS5LQUe1k6TrKuVyONVwWDWFe\nT9PUxzyLk0ZKhuIl9WIOFiGu+hH2tpwl5ZaabK4yeYB0ujywPF8TXzt7j57IVg0Fqsewm9dwBBVa\n7FUm0M0KirmumNhoyBkJq+/SL5iHymLuqvq73po6mE0EnjkMEYS8NP+g5wCOw+pJwCgVe9oVkxd+\nZqa5+yr8qZYlXOBrlvvyz8mbtNFMjQU5loaZxlC/zNWTf89qsmU5szvBifnyYK6oNspqi3257ddF\nF75apUSKXDEXDwJMXIEdbzlkRl+SUT8RRd9NrpZV7Irpe/LUnPjTyHIC9QRxnTLFru+1WtnmhTqm\n5My5hf4u48gcLo/wMi8/LfjN5bE3dxFFHIfZpuSDtNiPw+pJwBgVO5DSPHnWhDCG60GNidPk/GVX\nTJXOXXw9OYF6FCwvSWcO9av4s030kfpHXfVUR1CBL1K+hlmLvcL19LlNfxTeRODsSS/XN27eL08H\nhZEOfRX57Ri72djYVC6wWUckUA/mvAf4piZP80r3hlH1iBjNRKTdW/tXw3g9X21ll7tH5GuRy1O3\nS14hML1JP1iKwoo15q46v23PFZirjN7saabFPvRyAsAYFTulEzKKJsEcR1aZW1UDI/8eVLvK4SrM\niobmNaSPPblHHVGThTo4vlYTnAKLvcrDRlvqbfwoVln9DuUP6VdFg5gWe/wZV4gbe2+Y45XrsxN2\n+rhNjOCzlrCeHKzb510niWw6CiI8dcXH01eCuMQDUP696D5SpGiTxLLlffq71K4YXmFMJImIy+UD\nVnFi7oJo+bMQ1cvY3TTHbvJ0FXkJI3lWkkuJj72qHy6+B6WtOqB6eFwRegJVXj+5iKk46rpi9PlJ\nDelmPbIoKqZKR3eEjIWu6wJYh6L09GCVYo/buprf2NwfcWJR5vWn3ZlbKaqnbbKWsPZm1B05TIzI\nmKcPgjiM8mCRTGQ3sdgBFC5qrTfpZMNsgEEWLUfd3zUgHzxnMmG2gB45J++HzugUe17CSN734Igk\ni7HuUz378ND3Adb/0k3XQspiNy2uimFeJmbGY9PumOevrRqCqRV6H2sAa1dZtl7MSos9Z0Sy6qOZ\n5ywCLnRjuQ7BdaoV7moTR8hIooNFhOu2EjdaDdczgMS99PQVH37IuGF3gqeuBDhQi14DK9wjqk3K\nlG1eMhWQ5AYwSz/7qqgYvWsdxQ4gdwFvkfObBLTBFSHiCAz5Xm/LJrkR6s3lNWF0il26LtI+4LwO\nZw4J6/rYiTJLqKH+BGUWc9GEVCcysiHXiaElGA+hhn0qiQgyH5zVko4cgVqx000oqhejSwcXWY2J\nK4aTSKIV8qZcMRHH63gOiflU4OJBkIoCq2vM6O/PDxnbnoOtqYPtgHHpMIgzO8u+29hiL2mfopFW\nxNLwOfIj+MHqkat295XNw5iYylkr5KxyPgxCHPghwMDWwgVDAAQIEnDIgSMcuOTCFfJPbxMkIEio\nMFwBV/Sjcken2LMTm0U+YLNj1/U36phWk6blPPUEqh9yrismjsevbbEbk1vribYkp+nzr+pz3N1y\nK5eabUpRvZhVma+mK6YsJ8HEdMUEYYRZTrTNptmaClw8AA4XYbyt7lwHEcF15LKROmJq2xO4dCjr\nnAOrLHb5WuqKofz8g4gZUyGSDO2cUWJWOZ/cYQixwJXFsuVs/stgOMKJFfFUTHOV86EXYSZ8CCI8\n58QcW9NJrxZ4XUan2LOr4BTNoJuTeXWtFyr1sa//RU8cAT8Mi10xqD8iIDR3E5mYD7WSHKAlXEG1\nH6DrYia7mIQrIqDMPVUjkOLkmyACc3UrsU8mjlSKh4solm+d7+L6baku9O9l4og4YUxODhdfc9tz\n4K4IKxZgBEEAP6SU5Xzl6CrCqQxY2FeJfEQOJq58iJrKWSvl+aTcco7/r6icBUJccg7V555ANFyU\nvGtGp9jzXDF5nVhvW6eDa6vVHA3Erpg1ZNbICdTM5Kl6jeu91LyDlE8/FRoIZ1zPDL8EhjeZpK3R\nbOicnCgvm3TT3+XqyIv4HPW6KJk4HQJbUweXDgOp2NZ8wOY9tHY8B08FUSoCq8itIRzG/sI4OWMY\nLKIIixAQ5MUK2hUu9mc+Ts6ncIRMICIiXL/t4dS2V0s5N8G8fNNcjD4YnWI3XQWAdEOI3NoU8nUd\nxW5a0Uvxrg2+852ZzIY0rcqlB8carpj4/fqixZgr3TQN8ewKV8XNZ1Pgg4gxK1kFCEhGY1VdTPqY\nsjooQ2A+lW6Tw0W4tIhMljLlHHGUUqIhMw59X61Q5EMIkbKczb+UpZxjOT9z6OOpKwuc3d1OZaI+\n413B3szDxCUc+bJEt+dMe/NXA+1mT/fB6BS7OZlFSsnn5qIQYT51sFVh7cgs5lJZmjbqNAsibGd8\ntImPPf1/VczD21DA5sRx03mFLskuDqFjkFdZWzpzmeU/K+9jLqJRlFW5CfKUcxAdwQ8jhAzsLxap\nY82RYF3lfGrmAxA4sztr9PnjjHBmCCTvAYBEOou0b0+IHpUW1VEaGqNT7BptTZcVJzq9u8J8KyAd\nz51Y1F184bErJo6rXscVk75WE0wf+5DjeieC4oQhoHqpWgLiWjFVP5V2T5W5edahjuVcRTnfuDvH\nlSPgutkEZ3ZnjXzOJs+qsVRfGWapD12dO6m0KCdvdcntvnuc7jZ95GG0wegUe1zPBNWTTNa+Ryr6\nphvL1YzUMP+vfH7qWi1MnorEx940xLNLXEcmfOmRW1kBKhMixBZ71e9TH5bnhikNpcs48te1nLMK\nukg573khnrh4iL2tKban7SjjNskrKxBlEqqmrly7tP+69jI/YYh9PY/RKfa44TkpLdB2J0hb7Pp9\nOxbx8r0SXyOwXuap+doUHU5qppIPsbObi0NMHIozLvNW7zGRk+/pB7WZdJKnnA8DH0HImHkuLh9R\noXKeOJOUol5HOTdhNnFw3dYkrocyNFzDFaOJjTM1GJo6AocIN9LnRIOJ574ZnWLXRMZQuu1OkFeJ\nrqtynokrRt+7pitGv7Ykm75MyO2EeHaFnoAOQjms98MQYRQhinwcRVxoOR/4PkQgk5kmDmH/aLFS\nOVO0gB8CN5+cw5s4nSrnppze8TYtQiGOSBsxwLILclohGaorhLhGXDFE9McA3gBgAeAxAG9l5ott\nCLYuWukyTB9wu/fISz2vM3SvQxyTnRRVr0Vssbcmj/7s3Jmrq4hVlrP5IYMwwsEiwMVDByE7OPBl\n2qkjvKWYZ9NqnvBC+XGB+XRSaUJw3zvEVT/Cjtds8vBaJ2+RFNYjLdWu86mL3VlUqxxvW5ze8QY5\nOs2jqcV+P4B3MHNARH8E4B0Afru5WOuTWixBbWv76Z5Y7Kvj5ZuSDXesn3mqZ/PbkcccrawzeVpH\nOZf5nKfONE7RLnJreHQVO55Uzt+8dISQGWf38hdE18wm6nMxx9dZhUMya9gq9ebo9XE12YWlHUG4\nYXczo47ZZJgurDwaKXZm/pjx7wUAP91MnObon1bESSJx+66YvMnT6gs110WH05n3rnxu5rUpehEL\nPwixCAP4oY9DX0+mLitnAKlElKxyNq3nMp/zOpl+s0kYr7YTRFGlqBUCIQLHhaeqcP2Ot5SJbFkP\nWUs/7WM/DglBQ6NNH/vbAHyoxeuthWGwdxqOl12lp6vJUyAJp1vnY8TnZE7OWs4RR6nt2Q+jrWc/\nZCzCUJWoFZhNZtj15p0p5yZMHMLVhfxcQciYTVY3IBHAESvXWrUGl4rHKp82cATFyV6ASjC0I6Ha\nrFTsRPRxADfm7DrPzP+gjjkPIADw/pLr3APgHgC45ZZb1hK2CmZ4YJfheER5k6ft3wcongCtopwP\nfFmVbsqEy0dHsWtjXcs5ioAtcYDTux6u+iGO/AhntufdfPCGTB2B/ShAGFVfXCKJiuk/VtqSbzAd\nl0iUIbFSsTPza8r2E9GbAfw4gFdzyXiUme8DcB8AnDt3rrNxa1zmlrubPNXXTCv29UYG1ZSzH9cW\nv3y0qKScHeGAQLgyjTB3fMynLp59Yt7ccibldoly1mcdGK5yvRz6svpgpSE9DTs+f+y4KgEpzj/g\nequGWSRNo2LugJwsfQUzH7QjUjMSVwy3kuZffJ8kUSeIZChdyAGOgnBJOTPSoZDmpKCuPFemnEV0\nhCAEPNfFzafqKeco8jF1CFPHbaW2RrKGZnlW7xDQsey6XK1bYZUPYSQz2cnQ/jELuOnl+FbVtrEs\n0/SX/m4AHoD71Y/gAjPf21iqBqRcMVH12PJVlnNWOR/6PqIImHsLAAIhB3Bohpnr5SvnAtdGFbYm\nAkcIMXWd2so5qXFR67RSHLUggizbMFzlN1GK/GoNi908YsAfbbSY9WJcdFeqY+w0jYr5rrYEbacv\nQwAACyNJREFUaYOII0RRhCAKcDUIsQgiLIIIVxZBO5azoZxPej6eOQhx26kdBBHD5QOcuW6GHa/9\nnK8mkS1tx7EDyg8asbKmhvujE4LgCpEsMl2h8mLbtXUs9XDjJCX5f1ShcJtlmWOZeSp9byH2F/sp\n5ewKN7aGPceDS4S5K3Bme7sVy9lk1xO4fHiIq0FSt6Kr7qevv47bI44Db1E4PXFs1EAbLBOXECyk\ncqi0cIbxfshuprGi2zyIIkSRUGHEGxbqGHIsFft8MsfzTz6/UDlTeAW7MxdHamWXvVl5Uso6eK68\n9+EixLay0ruMYzdf1zq3RQ2sXTGMYfvYAe1XDytbfebHGfhHGyXmsoZJnRj7RdTlWCp2QQKeW5x9\npi3KiDn2s7YNEWE2ETj0w3gtz+7DHdc4V53T5m9DEMGPIhmKNnDtp1PPq0ycAukHYJsPQ0s19Mgq\nZF67VLUFGOV0M4FUqnq3VtfWxMEiiFSyTpcWe3NXTJsTUESIQ9KGbkzpOYDKFrvxi7D6ZDM4KlNY\nhysPeYJ+qIxTsRNk2d6Ow/F07YiDRQCgu2JYTSZARYNzi3COUUigttSr1vGxUTGbRwjpikn62IYF\nOoaMVrFLV0y3wzjPFRBEuOpH6r4dWezqdS2LvYNwR1OOoVvsU1dgdzbB3KtWwCkdFTPwDzdSXCEQ\nGtVDbVRMfUaq2NXkXseuAulnd+JEpa7u1SQqRhDgTRx4bnuV6Uw5hm6xA8ANu17lz5+OiulGHks5\nQgBhyEurJ1mqM0rFLgi9uQq2JnritLuyrUkhr3XOJZw9sRVP8LaBOQ85NuWXjooZ2Yc7Jjhq8tRa\n7OszSsVOSHzAXfeJ2VSv6NLdPZLJ0+7uUYe0K2YgQrVEajSyQTmuZVwh49f9KBpd/+qLUSp2WaCr\nn1Apz3XgCOr0Pkkc+zA6edoVs0FBOsbGT28GPSL0Q1uLfV1GqdhNU6uPJ/586sZrbHZBMnna2S1q\nkXbFDESolhjaQ/RaRCvzIIxGbTh0ybFMUFpFKsmkh45xemfa6fWbTJ52wZjdFXHc/4bluJaJC4FF\nHC9ebanHKBW76Nli79q6c9V6ml2sqboOY/axJ5m64/pcxwlnxP2rL0ap2E1FOwYfnec6uPX09qbF\niEk9OEfQviZdlDm21MP8zVrFvh6jHOfYWORukYXXup0w3hTj+0THDyKKlbv9/a7HKBV7YnV1F1t+\nrTNaxa5dMVajbBTdt8Yw4t4Eo1TsiP2kmxVjzBCN010Rl2DYsBzXOlqhW8NsPUap2O0EWPfI8qqb\nlqIb5Ehv01Jc2+hAAWuxr8coJ0/tBFj3zKcOVHLv6CBYo2DTCOtjb0QrFjsR/QYRMRGdbuN6TWlS\nDdFSjRPzKU5tdxu/vymIrCtm0+iQR/sbXo/Gip2IbgbwowC+1lycdrCuGEsTXEfAdUbppTw2OI5V\n7E1oo/e+C8BvARjMwFxnD9phnGUdnr03G+1o5LiwM3Vx/Y5nM0/XpFGrEdEbAfwvMz9U4dh7iOhB\nInrwySefbHLbCnLF9+z0PpZxYvvN5hGCsLc12bQYx5aVk6dE9HEAN+bsOg/gdwG8tsqNmPk+APcB\nwLlz5zq17rtYwNlisViOCysVOzO/Jm87Eb0YwK0AHlIWznMAfIaIXsrM/9eqlDVJXDFWs1sslmuP\ntcMdmfkRAGf0/0T0VQDnmPnbLcjVCGEnTy0WyzXMKGcmXEfg5HyK7YoLGFssFsuYaC1BiZmf19a1\n2uCkjWqwWCzXKKO02C0Wi+Vaxip2i8ViGRlWsVssFsvIsIrdYrFYRoZV7BaLxTIyrGK3WCyWkWEV\nu8VisYwMq9gtFotlZBBz/9V2iehJAP/d8DKnAWy8fEGGIcoEWLnqMESZACtXHYYoE9COXM9l5htW\nHbQRxd4GRPQgM5/btBwmQ5QJsHLVYYgyAVauOgxRJqBfuawrxmKxWEaGVewWi8UyMo6zYr9v0wLk\nMESZACtXHYYoE2DlqsMQZQJ6lOvY+tgtFovFks9xttgtFovFksPgFDsR3UFEXySiR4nod3L2e0T0\nIbX/U0T0PGPfO9T2LxLR63qW69eJ6D+J6GEi+hcieq6xLySiz6q/j/Qo01uI6Enj3j9v7HszEX1Z\n/b25LZkqyvUuQ6YvEdFFY19XbfVeIvoWEX2uYD8R0Z8qmR8mopcY+7psq1Vy/ayS52Ei+iQRfZ+x\n76tE9Ihqqwd7luuVRPSM8V39nrGv9PvvUKbfNOT5nOpLp9S+LtvqZiJ6gIi+QESfJ6JfzTmm3/7F\nzIP5A+AAeAzAbQCmAB4C8MLMMb8I4M/V+zcB+JB6/0J1vAe5FutjAJwe5XoVgLl6/wtaLvX//oba\n6i0A3p1z7ikAj6vXk+r9yb7kyhz/KwDe22Vbqev+MICXAPhcwf47AXwUAAH4AQCf6rqtKsr1Mn0/\nAK/Xcqn/vwrg9Iba65UA/qnp99+mTJlj3wDgEz211U0AXqLe7wL4Us5vsdf+NTSL/aUAHmXmx5l5\nAeCDAO7KHHMXgL9U7z8M4NVERGr7B5n5iJm/AuBRdb1e5GLmB5j5QP17AXJx7y6p0lZFvA7A/cz8\nFDM/DeB+AHdsSK6fAfCBlu5dCDP/G4CnSg65C8D7WHIBwAkiugndttVKuZj5k+q+QD/9qpJcJTTp\nl23K1Eu/AgBm/gYzf0a9vwzgCwDOZg7rtX8NTbGfBfA/xv9fx3IDxccwcwDgGQDXVzy3S7lM3g75\ndNbMiOhBIrpARHf3LNNPqaHfh4no5prndikXlLvqVgCfMDZ30VZVKJK7y7aqS7ZfMYCPEdF/ENE9\nG5DnB4noISL6KBG9SG3beHsR0RxSOf6NsbmXtiLpGr4dwKcyu3rtX62tedoSlLMtG7ZTdEyVc9el\n8rWJ6OcAnAPwCmPzLcz8BBHdBuATRPQIMz/Wg0z/COADzHxERPdCjnR+pOK5XcqleROADzNzaGzr\noq2qsIl+VRkiehWkYn+5sfmHVFudAXA/Ef2Xsmr74DOQ6e37RHQngL8H8AIMo73eAODfmdm07jtv\nKyLagXyY/BozX8ruzjmls/41NIv96wBuNv5/DoAnio4hIhfAHuTwrMq5XcoFInoNgPMA3sjMR3o7\nMz+hXh8H8K+QT/TOZWLm7xhy/AWA7696bpdyGbwJmeFyR21VhSK5u2yrShDR9wJ4D4C7mPk7ervR\nVt8C8Hdoz/W4Ema+xMz76v0/A5gQ0WkMoL1Q3q86aSsimkAq9fcz89/mHNJv/+piMqHBJIQLOXlw\nK5KJlxdljvklpCdP/0q9fxHSk6ePo73J0ypy3Q45afSCzPaTADz1/jSAL6OFyaSKMt1kvP8JABc4\nmbD5ipLtpHp/qq+2Usd9N+SEFnXdVsb1n4fiycAfQ3py69Ndt1VFuW6BnC96WWb7NoBd4/0nAdzR\no1w36u8OUkl+TbVdpe+/C5nUfm3obffVVupzvw/An5Qc02v/aq0TtNhId0LOKj8G4Lza9geQVjAA\nzAD8tersnwZwm3HueXXeFwG8vme5Pg7gmwA+q/4+ora/DMAjqoM/AuDtPcr0TgCfV/d+AMD3GOe+\nTbXhowDe2mdbqf9/H8AfZs7rsq0+AOAbAHxIK+ntAO4FcK/aTwD+TMn8CIBzPbXVKrneA+Bpo189\nqLbfptrpIfUdn+9Zrl82+tYFGA+evO+/D5nUMW+BDKIwz+u6rV4O6T552Pie7txk/7KZpxaLxTIy\nhuZjt1gsFktDrGK3WCyWkWEVu8VisYwMq9gtFotlZFjFbrFYLCPDKnaLxWIZGVaxWywWy8iwit1i\nsVhGxv8DCmw+u2txfdEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111689a90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as st\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# Generate data\n",
    "M = 20\n",
    "\n",
    "X = np.linspace(0, 2, num=M)\n",
    "# X = np.array([5, 14, 19], dtype=np.float)\n",
    "t_real = np.sin(X)\n",
    "t = t_real + np.random.randn(M) * 0.25\n",
    "\n",
    "plt.scatter(X, t, label='Data points')\n",
    "\n",
    "# Generative process: p(t|W,X,beta) = N(t|XW+b,beta)\n",
    "beta = 1\n",
    "\n",
    "# Prior: N(w|0,1/alpha*I)\n",
    "alpha = 1\n",
    "\n",
    "# Posterior: N(w|m,s):\n",
    "s = 1/(alpha + beta * X.T @ X)\n",
    "m = beta * s * X.T @ t\n",
    "\n",
    "# Infer p(t|t,alpha,beta) the predictive distribution\n",
    "X_pred = np.linspace(0, 2, num=100)\n",
    "\n",
    "m_pred = m * X_pred\n",
    "s_pred = 1/beta + X_pred.T @ X_pred * s\n",
    "std_pred = np.sqrt(s_pred)\n",
    "\n",
    "plt.plot(X_pred, m_pred, color='red', alpha=0.75, label='Regression line')\n",
    "plt.fill_between(\n",
    "    X_pred, m_pred-std_pred, m_pred+std_pred,\n",
    "    interpolate=True, color='green', alpha=0.1, label='+- 1 stddev'\n",
    ")\n",
    "\n",
    "# Sample from predictive dist.\n",
    "ys = np.random.normal(m_pred, std_pred)\n",
    "\n",
    "plt.plot(X_pred, ys, alpha=0.15, label='Posterior samples')\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# #2 Bayesian Dark Knowledge\n",
    "\n",
    "Bayesian Dark Knowledge\n",
    "Korattikara, et. al., 2015"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Successfully downloaded train-images-idx3-ubyte.gz 9912422 bytes.\n",
      "Extracting ../../MNIST_data/train-images-idx3-ubyte.gz\n",
      "Successfully downloaded train-labels-idx1-ubyte.gz 28881 bytes.\n",
      "Extracting ../../MNIST_data/train-labels-idx1-ubyte.gz\n",
      "Successfully downloaded t10k-images-idx3-ubyte.gz 1648877 bytes.\n",
      "Extracting ../../MNIST_data/t10k-images-idx3-ubyte.gz\n",
      "Successfully downloaded t10k-labels-idx1-ubyte.gz 4542 bytes.\n",
      "Extracting ../../MNIST_data/t10k-labels-idx1-ubyte.gz\n",
      "Iter: 1000; S_loss: 690.7835\n",
      "Iter: 2000; S_loss: 690.7092\n",
      "Iter: 3000; S_loss: 690.7138\n",
      "Iter: 4000; S_loss: 690.7006\n",
      "Iter: 5000; S_loss: 690.7684\n",
      "Iter: 6000; S_loss: 690.7776\n",
      "Iter: 7000; S_loss: 690.7752\n",
      "Iter: 8000; S_loss: 690.7309\n",
      "Iter: 9000; S_loss: 690.8474\n",
      "Iter: 10000; S_loss: 690.6991\n",
      "Test accuracy: 0.0464\n"
     ]
    }
   ],
   "source": [
    "\n",
    "import numpy as np\n",
    "from tensorflow.examples.tutorials.mnist import input_data\n",
    "\n",
    "\n",
    "mnist = input_data.read_data_sets('../../MNIST_data', one_hot=True)\n",
    "X_dim = mnist.train.images.shape[1]\n",
    "y_dim = mnist.train.labels.shape[1]\n",
    "\n",
    "N = mnist.train.images.shape[0]  # training set size\n",
    "M = 100  # minibatch size\n",
    "T = 10000  # num of iteration\n",
    "\n",
    "Theta = np.vstack([np.zeros([1, y_dim]), np.random.randn(X_dim, y_dim) * 0.001])\n",
    "W = np.vstack([np.zeros([1, y_dim]), np.random.randn(X_dim, y_dim) * 0.001])\n",
    "\n",
    "\n",
    "def softmax(x):\n",
    "    m = np.max(x)\n",
    "    return np.exp(x - m) / np.sum(np.exp(x - m))\n",
    "\n",
    "\n",
    "# SGLD Teacher-Student distillation\n",
    "# ---------------------------------\n",
    "\n",
    "eta = 0.1  # teacher learning rate\n",
    "rho = 0.1  # student learning rate\n",
    "lam = 1  # teacher prior\n",
    "gam = 1e-3  # student prior\n",
    "\n",
    "burnin = 1000\n",
    "thinning = 10\n",
    "\n",
    "for t in range(1, T+1):\n",
    "    # Train teacher\n",
    "    # -------------\n",
    "    X_mb, y_mb = mnist.train.next_batch(M)\n",
    "    X_train = np.hstack([X_mb, np.ones([M, 1])])  # add bias dim\n",
    "\n",
    "    eta_t = eta/t\n",
    "\n",
    "    p = softmax(X_train @ Theta)\n",
    "    grad_p = p - y_mb\n",
    "\n",
    "    grad_loglik = X_train.T @ grad_p  # 784x16 . 16x10\n",
    "    grad_logprior = lam * Theta\n",
    "    grad_logpost = grad_logprior + N/M * grad_loglik\n",
    "\n",
    "    z = np.random.normal(0, np.sqrt(eta_t))\n",
    "    delta = eta_t/2 * grad_logpost + z\n",
    "    Theta += delta\n",
    "\n",
    "    if t > burnin and t % thinning == 0:\n",
    "        # Train student\n",
    "        # -------------\n",
    "        X_s_mb = X_mb + np.random.normal(0, 1e-3)  # perturb\n",
    "        X_s_train = np.hstack([X_s_mb, np.ones([M, 1])])  # add bias dim\n",
    "        rho_t = rho/t\n",
    "\n",
    "        s = softmax(X_s_train @ W)\n",
    "        p = softmax(X_s_train @ Theta)\n",
    "        grad_s = s - p\n",
    "\n",
    "        grad_loglik = X_s_train.T @ grad_s  # 784x16 . 16x10\n",
    "        grad_logprior = gam * W\n",
    "        grad_logpost = grad_logprior + 1/M*grad_loglik\n",
    "\n",
    "        delta = rho_t * grad_logpost\n",
    "        W -= delta\n",
    "\n",
    "    # Diagnostics\n",
    "    # -----------\n",
    "    if t % 1000 == 0:\n",
    "        s = softmax(X_train @ W)\n",
    "        loss = -gam/2 * np.sum(W**2) - np.sum(y_mb * np.log(s + 1e-8))\n",
    "        print('Iter: {}; S_loss: {:.4f}'.format(t, loss))\n",
    "\n",
    "\n",
    "# Test\n",
    "# ----\n",
    "\n",
    "X_test, y_test = mnist.test.images, mnist.test.labels\n",
    "X_test = np.hstack([X_test, np.ones([X_test.shape[0], 1])])  # add bias dim\n",
    "\n",
    "y = softmax(X_test @ W)\n",
    "acc = np.mean(y.argmax(axis=1) == y_test.argmax(axis=1))\n",
    "\n",
    "print('Test accuracy: {:.4f}'.format(acc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# #3 GMM Gibbs \n",
    "\n",
    "Posterior sampling for Gaussian Mixture Model using Gibbs sampler\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20 data in cluster-0, mean: [ 4.51710905  5.00999685]\n",
      "20 data in cluster-1, mean: [ 7.77123298  7.78397887]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as st\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "K = 2\n",
    "\n",
    "# Generate data\n",
    "X1 = np.random.multivariate_normal([5, 5], np.diag([0.5, 0.5]), size=20)\n",
    "X2 = np.random.multivariate_normal([8, 8], np.diag([0.5, 0.5]), size=20)\n",
    "X = np.vstack([X1, X2])\n",
    "\n",
    "N = X.shape[0]\n",
    "\n",
    "# GMM params\n",
    "mus = np.array([[1, 1], [15, 15]], dtype='float')\n",
    "sigmas = np.array([np.diag([1, 1]), np.diag([1, 1])], dtype='float')\n",
    "lambdas = np.array([np.linalg.inv(sigmas[0]), np.linalg.inv(sigmas[1])])\n",
    "pis = np.array([0.5, 0.5])  # Mixing probs.\n",
    "zs = np.zeros([N])  # Assignments\n",
    "\n",
    "# Priors\n",
    "alpha = np.array([1, 1])\n",
    "pis = np.random.dirichlet(alpha)\n",
    "mus0 = np.array([[1, 1], [1, 1]], dtype='float')\n",
    "sigmas0 = np.array([np.diag([1, 1]), np.diag([1, 1])], dtype='float')\n",
    "lambdas0 = np.array([np.linalg.inv(sigmas0[0]), np.linalg.inv(sigmas0[1])])\n",
    "\n",
    "# Gibbs sampler\n",
    "for it in range(50):\n",
    "    # Sample from full conditional of assignment\n",
    "    # z ~ p(z) \\propto pi*N(y|pi)\n",
    "    probs = np.zeros([N, K])\n",
    "\n",
    "    for k in range(K):\n",
    "        p = pis[k] * st.multivariate_normal.pdf(X, mean=mus[k], cov=sigmas[k])\n",
    "        probs[:, k] = p\n",
    "\n",
    "    # Normalize\n",
    "    probs /= np.sum(probs, axis=1)[:, np.newaxis]\n",
    "\n",
    "    # For each data point, draw the cluster assignment\n",
    "    for i in range(N):\n",
    "        z = np.random.multinomial(n=1, pvals=probs[i]).argmax()\n",
    "        zs[i] = z\n",
    "\n",
    "    # Sample from full conditional of cluster parameter\n",
    "    # Assume fixed covariance => posterior is Normal\n",
    "    # mu ~ N(mu, sigma)\n",
    "    Ns = np.zeros(K, dtype='int')\n",
    "\n",
    "    for k in range(K):\n",
    "        # Gather all data points assigned to cluster k\n",
    "        Xk = X[zs == k]\n",
    "        Ns[k] = Xk.shape[0]\n",
    "\n",
    "        # Covariance of posterior\n",
    "        lambda_post = lambdas0[k] + Ns[k]*lambdas[k]\n",
    "        cov_post = np.linalg.inv(lambda_post)\n",
    "\n",
    "        # Mean of posterior\n",
    "        left = cov_post\n",
    "        right = lambdas0[k] @ mus0[k] + Ns[k]*lambdas[k] @ np.mean(Xk, axis=0)\n",
    "        mus_post = left @ right\n",
    "\n",
    "        # Draw new mean sample from posterior\n",
    "        mus[k] = st.multivariate_normal.rvs(mus_post, cov_post)\n",
    "\n",
    "    # Sample from full conditional of the mixing weight\n",
    "    # pi ~ Dir(alpha + n)\n",
    "    pis = np.random.dirichlet(alpha + Ns)\n",
    "\n",
    "# Expected output:\n",
    "# ----------------\n",
    "# 20 data in cluster-0, mean: [ 5  5 ]\n",
    "# 20 data in cluster-1, mean: [ 8  8 ]\n",
    "for k in range(K):\n",
    "    print('{} data in cluster-{}, mean: {}'.format(Ns[k], k, mus[k]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# #4 GMM Variation\n",
    "\n",
    "Posterior sampling for 2D-Gaussian Mixture Model using Mean-Field VI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20 data in cluster-0, mean: [-4.45036297 -4.87464662]\n",
      "20 data in cluster-1, mean: [ 17.07966899  17.1744302 ]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as st\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "K = 2\n",
    "\n",
    "# Generate data\n",
    "X1 = np.random.multivariate_normal([-5, -5], np.diag([2, 0.5]), size=20)\n",
    "X2 = np.random.multivariate_normal([18, 18], np.diag([0.7, 1]), size=20)\n",
    "X = np.vstack([X1, X2])\n",
    "\n",
    "N = X.shape[0]\n",
    "\n",
    "# GMM params initialization\n",
    "mu = np.array([[1, 1], [2, 2]], dtype=float)\n",
    "var = np.array([1, 1], dtype=float)\n",
    "phi = np.zeros([N, K]) + 1/K\n",
    "c = np.random.randint(K, size=N)  # Assignments\n",
    "\n",
    "# Priors\n",
    "var0 = 1\n",
    "\n",
    "for it in range(5):\n",
    "    # Update variational param phi, the assignment probs\n",
    "    for k in range(K):\n",
    "        phi[:, k] = np.exp(X @ mu[k] - (2*var[k] + mu[k].T @ mu[k])/2)\n",
    "\n",
    "    # Normalize\n",
    "    phi /= np.sum(phi, axis=1)[:, np.newaxis]\n",
    "    # Update assignments\n",
    "    c = np.argmax(phi, axis=1)\n",
    "\n",
    "    # Update variational param mu and var, the params of Gaussian component\n",
    "    for k in range(K):\n",
    "        sum_phi = np.sum(phi[:, k])\n",
    "        mu[k] = phi[:, k] @ X / (1/var0 + sum_phi)\n",
    "        var[k] = 1 / (1/var0 + sum_phi)\n",
    "\n",
    "\n",
    "# Expected output:\n",
    "# ----------------\n",
    "# 20 data in cluster-0, mean: [ -5  -5 ]\n",
    "# 20 data in cluster-1, mean: [ 18  18 ]\n",
    "for k in range(K):\n",
    "    n = np.sum(c == k)\n",
    "    print('{} data in cluster-{}, mean: {}'.format(n, k, mu[k]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# #5 GMM Variational 1-D \n",
    "\n",
    "Posterior sampling for Gaussian Mixture Model using Gibbs sampler\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "76 data in cluster-0, mean: 2.8807925918438864\n",
      "4 data in cluster-1, mean: -4.795346463816408\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as st\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "K = 2\n",
    "\n",
    "# Generate data\n",
    "X1 = np.random.normal(-5, 0.5, size=20)\n",
    "X2 = np.random.normal(3, 0.5, size=20)\n",
    "X = np.concatenate([X1, X2])\n",
    "\n",
    "N = X.shape[0]\n",
    "\n",
    "# GMM params initialization\n",
    "mu = np.array([50, -20], dtype=float)\n",
    "var = np.array([1, 1], dtype=float)\n",
    "phi = np.zeros([N, K]) + 1/K\n",
    "c = np.zeros([N, K])  # Assignments\n",
    "\n",
    "# Priors\n",
    "mu0 = 0\n",
    "var0 = 1\n",
    "\n",
    "for it in range(5):\n",
    "    # Update variational param phi, the assignment probs\n",
    "    for k in range(K):\n",
    "        phi[:, k] = np.exp(mu[k]*X - (mu[k]**2 + var[k])/2)\n",
    "\n",
    "    # Normalize\n",
    "    phi /= np.sum(phi, axis=1)[:, np.newaxis]\n",
    "    # Update assignments\n",
    "    c[np.argmax(phi, axis=1)] = 1\n",
    "\n",
    "    # Update variational param mu and var, the params of Gaussian component\n",
    "    for k in range(K):\n",
    "        sum_phi = np.sum(phi[:, k])\n",
    "        mu[k] = phi[:, k] @ X / (1/var0 + sum_phi)\n",
    "        var[k] = 1 / (1/var0 + sum_phi)\n",
    "\n",
    "\n",
    "# Expected output:\n",
    "# ----------------\n",
    "# 20 data in cluster-0, mean: -5\n",
    "# 20 data in cluster-1, mean: 3\n",
    "for k in range(K):\n",
    "    n = np.sum(c == k)\n",
    "    print('{} data in cluster-{}, mean: {}'.format(n, k, mu[k]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# #6 LDA with Collapsed Gibb's Sampler \n",
    "\n",
    "LDA with Collapsed Gibbs Sampler\n",
    "\n",
    "Reference: Kevin Murphy's book Ch. 27\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Documents:\n",
      "----------\n",
      "[[0 0 1 2 2]\n",
      " [0 0 1 1 1]\n",
      " [0 1 2 2 2]\n",
      " [4 4 4 4 4]\n",
      " [3 3 4 4 4]\n",
      " [3 4 4 4 4]]\n",
      "\n",
      "Document topic distribution:\n",
      "----------------------------\n",
      "[[ 0.14285714  0.85714286]\n",
      " [ 0.57142857  0.42857143]\n",
      " [ 0.57142857  0.42857143]\n",
      " [ 0.71428571  0.28571429]\n",
      " [ 0.42857143  0.57142857]\n",
      " [ 0.28571429  0.71428571]]\n",
      "\n",
      "Topic word distribution:\n",
      "----------------------------\n",
      "[[ 0.33333333  0.22222222  0.11111111  0.05555556  0.27777778]\n",
      " [ 0.09090909  0.18181818  0.27272727  0.31818182  0.13636364]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "# Words\n",
    "W = np.array([0, 1, 2, 3, 4])\n",
    "\n",
    "# X := document words\n",
    "X = np.array([\n",
    "    [0, 0, 1, 2, 2],\n",
    "    [0, 0, 1, 1, 1],\n",
    "    [0, 1, 2, 2, 2],\n",
    "    [4, 4, 4, 4, 4],\n",
    "    [3, 3, 4, 4, 4],\n",
    "    [3, 4, 4, 4, 4]\n",
    "])\n",
    "\n",
    "N_D = X.shape[0]  # num of docs\n",
    "N_W = W.shape[0]  # num of words\n",
    "N_K = 2  # num of topics\n",
    "\n",
    "# Dirichlet priors\n",
    "alpha = 1\n",
    "gamma = 1\n",
    "\n",
    "\n",
    "# --------------\n",
    "# Initialization\n",
    "# --------------\n",
    "\n",
    "# Z := word topic assignment\n",
    "Z = np.zeros(shape=[N_D, N_W], dtype=int)\n",
    "\n",
    "# Some counts for sufficient statistics\n",
    "c_ik = np.zeros(shape=[N_D, N_K], dtype=int)\n",
    "c_vk = np.zeros(shape=[N_W, N_K], dtype=int)\n",
    "c_k = np.zeros(shape=N_K, dtype=int)\n",
    "\n",
    "for i in range(N_D):\n",
    "    for v in range(N_W):\n",
    "        # Randomly assign word's topic\n",
    "        k = np.random.randint(N_K)\n",
    "        Z[i, v] = k\n",
    "\n",
    "        # Record counts\n",
    "        c_ik[i, k] += 1\n",
    "        c_vk[v, k] += 1\n",
    "        c_k[k] += 1\n",
    "\n",
    "\n",
    "L = np.array([x.size for x in X])\n",
    "\n",
    "\n",
    "# --------------\n",
    "# Gibbs sampling\n",
    "# --------------\n",
    "\n",
    "for it in range(1000):\n",
    "    # Sample from full conditional of Z\n",
    "    # ---------------------------------\n",
    "    for i in range(N_D):\n",
    "        for v in range(N_W):\n",
    "            # Sufficient statistics for the full conditional\n",
    "            k = Z[i, v]\n",
    "            c_ik[i, k] -= 1\n",
    "            c_vk[v, k] -= 1\n",
    "            c_k[k] -= 1\n",
    "\n",
    "            # Calculate full conditional p(z_iv | .)\n",
    "            left = (c_vk[v, :] + gamma) / (c_k + N_W*gamma)\n",
    "            right = (c_ik[i, :] + alpha) / (L[i] + N_K*alpha)\n",
    "\n",
    "            p_z_iv = left * right\n",
    "            p_z_iv /= np.sum(p_z_iv)\n",
    "\n",
    "            # Resample word topic assignment\n",
    "            k = np.random.multinomial(1, p_z_iv).argmax()\n",
    "\n",
    "            # Update counts\n",
    "            Z[i, v] = k\n",
    "            c_ik[i, k] += 1\n",
    "            c_vk[v, k] += 1\n",
    "            c_k[k] += 1\n",
    "\n",
    "\n",
    "print('Documents:')\n",
    "print('----------')\n",
    "print(X)\n",
    "\n",
    "print()\n",
    "\n",
    "print('Document topic distribution:')\n",
    "print('----------------------------')\n",
    "print((c_ik + alpha) / np.sum(c_ik + alpha, axis=1)[:, np.newaxis])\n",
    "\n",
    "print()\n",
    "\n",
    "print('Topic word distribution:')\n",
    "print('----------------------------')\n",
    "x = (c_vk + gamma).T\n",
    "print(x / np.sum(x, axis=1)[:, np.newaxis])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# #7 LDA with Gibb's Sampler \n",
    "\n",
    "LDA with Gibbs Sampler\n",
    "Reference: Kevin Murphy's book Ch. 27\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Documents:\n",
      "----------\n",
      "[[0 0 1 2 2]\n",
      " [0 0 1 1 1]\n",
      " [0 1 2 2 2]\n",
      " [4 4 4 4 4]\n",
      " [3 3 4 4 4]\n",
      " [3 4 4 4 4]]\n",
      "\n",
      "Document topic distribution:\n",
      "----------------------------\n",
      "[[ 0.53798222  0.46201778]\n",
      " [ 0.86426942  0.13573058]\n",
      " [ 0.78028235  0.21971765]\n",
      " [ 0.02622966  0.97377034]\n",
      " [ 0.0618754   0.9381246 ]\n",
      " [ 0.22269837  0.77730163]]\n",
      "\n",
      "Topic's word distribution:\n",
      "-------------------------\n",
      "[[ 0.21615902  0.23579512  0.45697282  0.08111374  0.0099593 ]\n",
      " [ 0.02411378  0.02508557  0.05772617  0.19206773  0.70100675]]\n",
      "\n",
      "Word topic assignment:\n",
      "-------------------------\n",
      "[[ 0.  0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.]\n",
      " [ 1.  1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.  1.]\n",
      " [ 0.  1.  1.  1.  1.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "# Words\n",
    "W = np.array([0, 1, 2, 3, 4])\n",
    "\n",
    "# D := document words\n",
    "X = np.array([\n",
    "    [0, 0, 1, 2, 2],\n",
    "    [0, 0, 1, 1, 1],\n",
    "    [0, 1, 2, 2, 2],\n",
    "    [4, 4, 4, 4, 4],\n",
    "    [3, 3, 4, 4, 4],\n",
    "    [3, 4, 4, 4, 4]\n",
    "])\n",
    "\n",
    "N_D = X.shape[0]  # num of docs\n",
    "N_W = W.shape[0]  # num of words\n",
    "N_K = 2  # num of topics\n",
    "\n",
    "# Dirichlet priors\n",
    "alpha = 1\n",
    "gamma = 1\n",
    "\n",
    "\n",
    "# --------------\n",
    "# Initialization\n",
    "# --------------\n",
    "\n",
    "# Z := word topic assignment\n",
    "Z = np.zeros(shape=[N_D, N_W])\n",
    "for i in range(N_D):\n",
    "    for l in range(N_W):\n",
    "        Z[i, l] = np.random.randint(N_K)  # randomly assign word's topic\n",
    "\n",
    "# Pi := document topic distribution\n",
    "Pi = np.zeros([N_D, N_K])\n",
    "for i in range(N_D):\n",
    "    Pi[i] = np.random.dirichlet(alpha*np.ones(N_K))\n",
    "\n",
    "# B := word topic distribution\n",
    "B = np.zeros([N_K, N_W])\n",
    "for k in range(N_K):\n",
    "    B[k] = np.random.dirichlet(gamma*np.ones(N_W))\n",
    "\n",
    "\n",
    "# --------------\n",
    "# Gibbs sampling\n",
    "# --------------\n",
    "\n",
    "for it in range(1000):\n",
    "    # Sample from full conditional of Z\n",
    "    # ---------------------------------\n",
    "    for i in range(N_D):\n",
    "        for l in range(N_W):\n",
    "            # Calculate params for Z\n",
    "            p_bar_il = np.exp(np.log(Pi[i]) + np.log(B[:, X[i, l]]))\n",
    "            p_il = p_bar_il / np.sum(p_bar_il)\n",
    "\n",
    "            # Resample word topic assignment Z\n",
    "            z_il = np.random.multinomial(1, p_il)\n",
    "            Z[i, l] = np.argmax(z_il)\n",
    "\n",
    "    # Sample from full conditional of Pi\n",
    "    # ----------------------------------\n",
    "    for i in range(N_D):\n",
    "        m = np.zeros(N_K)\n",
    "\n",
    "        # Gather sufficient statistics\n",
    "        for k in range(N_K):\n",
    "            m[k] = np.sum(Z[i] == k)\n",
    "\n",
    "        # Resample doc topic dist.\n",
    "        Pi[i, :] = np.random.dirichlet(alpha + m)\n",
    "\n",
    "    # Sample from full conditional of B\n",
    "    # ---------------------------------\n",
    "    for k in range(N_K):\n",
    "        n = np.zeros(N_W)\n",
    "\n",
    "        # Gather sufficient statistics\n",
    "        for v in range(N_W):\n",
    "            for i in range(N_D):\n",
    "                for l in range(N_W):\n",
    "                    n[v] += (X[i, l] == v) and (Z[i, l] == k)\n",
    "\n",
    "        # Resample word topic dist.\n",
    "        B[k, :] = np.random.dirichlet(gamma + n)\n",
    "\n",
    "\n",
    "print('Documents:')\n",
    "print('----------')\n",
    "print(X)\n",
    "\n",
    "print()\n",
    "\n",
    "print('Document topic distribution:')\n",
    "print('----------------------------')\n",
    "print(Pi)\n",
    "\n",
    "print()\n",
    "\n",
    "print('Topic\\'s word distribution:')\n",
    "print('-------------------------')\n",
    "print(B)\n",
    "\n",
    "print()\n",
    "\n",
    "print('Word topic assignment:')\n",
    "print('-------------------------')\n",
    "print(Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}