{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Variational Coin Toss\n", "\n", "Code to go along with [this blog post](http://www.openias.org/variational-coin-toss)." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtwAAAGuCAYAAACjq9E+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xls3PWd//Hn2B6fie/YcWzHzg1OHBITCDkcjiSQ7bIL\nS4GURoIWUYUKurC0rKi21ZZtUVoBisoWsVDUCtAmheZHN0s3pA3hSCAHpAnEuZzLju3Y8RXf19ie\n+f3xxfZ8Z8Y5PfOd4/WQrHTe/o7zdgHP2995f95vm8vlciEiIiIiIn4RZXUCIiIiIiLhTAW3iIiI\niIgfqeAWEREREfEjFdwiIiIiIn6kgltERERExI9irE7An3p7ezl06BATJkwgOjra6nREREREJAwN\nDg7S2NjInDlziI+P9/p8WBfchw4dYs2aNVanISIiIiIR4L//+79ZsGCBVzysC+4JEyYAxjc/ceJE\ni7MRERERkXB07tw51qxZM1x7egrrgnuojWTixInk5eVZnI2IiIiIhLPRWph1aFJERERExI9UcIuI\niIiI+JEKbhERERERP1LBLSIiIiLiRyq4RURERET8SAW3iIiIiIgfqeAWEREREfEjFdwiIiIiIn6k\ngltERERExI9UcIuIiIiI+JEKbhERERERP1LBLSIiIiLiRzFWJyAiIv7XPwgdfdA+9NELnQ7oGzQ+\n5xgAxyD0OyHKBjYbRNuM/x0dBQl2SHT7SIqFtARITwB7tNXfnYhIcFPBLSISZhwD0NBlfDR2Gn+2\n9vrv70uOM4rvjETIGW98JNj99/eJiIQaFdwiIiFu0AnnOqCqDapajQI7kIbump9pHYmlJ8CkZMhN\nhsmpEK9XGxGJYPoRKCISgvoG4PR5OHUeatqMVpBgcr7H+DhUDzaMwntqOkxJN+6Ii4hEEhXcIiIh\non8QKlvgeBNUtoLTdXnPHx9nFLvJcZAcb/wZHwOx0UYf9tCfLpfxtYc++gehux96+o0/u/uhrRda\neow/L5aGC6hpNz52VMKEJLh2AszMVOuJiEQGFdwiIkGuqQvK6qG88dLvZKclQFYSZI0z/sxMMgrq\nsTbghNYeaO6Gc51Q2w5N3Rd+TmOX8fHpGZiSBtdmQUGqcUBTRCQcqeAWEQlCg06jXeTgOajruPj1\niXajV3pyCuSnQGKs/3MEiIkyivnMJJg1wYj1DRg95TXtxh358z2+n+t0Gd/jqfPG1JPrJsKcbIjT\nK5OIhBn9WBMRCSKOQSg7B1/WGa0bF5IcZ7RlTM+AzERjlF8wiIuBgjTjY0mBcQf8dAtUnDd+efDV\ngtLlgF1V8EUNzM6G63LU6y0i4UMFt4hIEOgbgK/OwVd10Dsw+nVJdpiRaRTaWUnBU2RfSGoClCRA\nySRj9vexRjja4HtUYb/T+GXjqzrj+7wxz2iPEREJZSq4RUQs1DcAB2qNYtsxOPp1eclQPNHoeY4O\n4R3B42JhQS5cP8m42320AY43G73g7lwYh0NPNBk93jfmGYc+RURCkQpuERELDDrhcAPsrR79jrY9\nyig2i7MhPTGw+fmbzWbM6Z6UDIsLjPGBB895t9G4gCMNxl3x4myjWA9Uf7qIyFhRwS0iEkAuF1S0\nwGdnRt/+GBtt9DDPmwjxETA2L8EON+TB/EnGJJYDdcbIQXdOl/EuwJFG4273dRND+06/iEQWFdwi\nIgFyvhs+qTCmd/gSHwPzcmDuxMic1BETZRyYLMqCE83G3X/PX0r6B41fVo40wLJCYzKLiEiwi8Af\n6SIigTXghH018Lda38tqYqKMA4XzJ/lnVnaosdlGpq8cbYDPa4zDlu5aemDzUZiWDksLNdFERIKb\nCm4RET+qboOPT4/ePlKUBTflG3OoxSzKZtzxnjXB6PHeW+19sPTUeahqNfrAi7NDY2qLiEQeFdwi\nIn7QN2CsMT/W6Pvz+SmwtMBYGCMXFhNltNrMzITdVUY7ibt+p9Gqc7IZlk+DlHhr8hQRGY0KbhGR\nMVbdBh+c9G6DAGMj5LIpMD1dd2MvV6LdKKhnZxkFdkOX+fNn22HDV7BosnGoUv//ikiwUMEtIjJG\nBpzGHdgv63x/vjjbKAYj8UDkWJo4Hu4vNsYqfnrGOEg5ZMAJOyvhVDPcPkOzu0UkOOjHvojIGGjs\ngr+egPM93p9LT4DbpkHO+MDnFa5sNpiTDQWp8OEpqGozf762AzZ+ZdwRn5ZhTY4iIkNUcIuIXAWX\nC8rqjbuqviaQlEwyDkVqZrR/jI+Df7wWjjYa/wzcD1X2DcKW48Y7C0sLIEYTYETEIiq4RUSukGPQ\nuLt6otn7c+PjYOV0yE0OfF6RxmYzpr1MToHtPu52l9VDbTusmhl+GztFJDQEtOAuLy/nhz/8Id3d\n3Xz44Yc+r/nJT37C5s2bTTGn00lJSQlvvfUWzzzzDJs3byYmZiT1mJgYDhw44NfcRUTcNXXB+8d9\nj/u7doKxlCVWtzQCatzXd7sP1Bm99O7vODT3wDtlsGK6Md9bRCSQAvZysGXLFtatW8fcuXM5evTo\nqNf94he/4Be/+MXwY5fLxQMPPMCdd945HLvrrrv45S9/6dd8RURGc7TRmK094DTH7dGwfCrMyLQm\nLzHudpdMMt5Z+MtxaOsb+Vy/0/glqWSScXg1SlNMRCRAAtZV2N3dzdtvv82iRYsu63mbNm2iv7+f\n++67z0+ZiYhcGqfL6BP+4KR3sZ2RCKuLVWwHi+xx8K25MMvHP4/9tfC/R6GnP/B5iUhkCljBfe+9\n9zJp0qTLek5PTw/r16/nxz/+MVFRI6mWl5ezevVqSkpKuPvuu9m/f/9YpysiYtI3AO8d8z3yr2gC\n3D8H0hICn5eMLjbG6KO/dar33ezqNqPFpLHL93NFRMZSUJ+b37hxI9OmTWPBggXDsfz8fAoLC3nx\nxRfZuXMnpaWlPPLIIzQ3+zi1JCIyBlp74I+HjBXi7qJtsGIaLJ+uCRjBamh84D2zIclu/lx7H/y/\nQ1DRYk1uIhI5grbgdjqd/P73v+fhhx82xR977DHWr19PXl4eSUlJPPnkkyQmJrJt2zaLMhWRcDZ0\nJ7TFY752UizcOweuzbImL7k8OeNh9VzvWej9Tvi/Y3DwnDV5iUhkCNqCe9++fXR2drJkyZILXhcd\nHU1OTg4NDQ0BykxEIsWxRqPXt2/QHM8eZ/RrZ42zJi+5Mkmx8E9Fxlxudy6MVfGjzVIXEblaQVtw\nf/DBB9x0003ExsYOx1wuF+vWrePYsWPDMYfDQVVVFfn5+VakKSJhyOWCv52FbSe9C7BZmV+3J8T6\nfq4Et+gouGWqMbbRc0jJl3Wwpdy8Kl5EZCwERcG9atUq9u7da4odOXKEvLw8U8xms1FTU8Ozzz5L\nfX09XV1dvPDCC9jtdm6//fZApiwiYcrpgh2VsKvK+3OLJhuH8GKC4ienXI3rcuDvrwG7xz/Lihb4\n0xFNMBGRsRWwOdx33HEHtbW1OJ1OBgYGKC4uBmDr1q1UVFTQ3d1tur6xsZFly5Z5fZ3nnnuOX/3q\nV9xzzz10dnYyd+5c3njjDZKSkgLyfYhI+BpwwrYTcPK8OR5lMwrtmRr5F1ampME358B7R6HLrcCu\n74R3D8Nd1xrLdERErpbN5XKFbcdaTU0Ny5cvZ/v27V53y0VE3DkG4c/H4Gy7OR4bDX8/C/JSrMlL\n/K+zzxj52GS+78P4OLj7WkjVuEcRuYiL1Zx6Y1REIl7vAPzPEe9iO8kO35ytYjvcjYsz/jnnJpvj\nHX2w6bBmdYvI1VPBLSIRrbsf/nTYaCNwl5YA9xVDprrVIkJsDPzjtUabibuefqO9pLbd9/NERC6F\nCm4RiVidfUYx5dlKkD0O7p1ttBRI5IiJgm/MgmsnmOOOQWM85Nk2a/ISkdCngltEIlJ7L/y/w94L\nbSaNN/p24+2+nyfhLcoGy6fBvBxzvN8J/3sMalR0i8gVUMEtIhGnvde4s93eZ45PTjHaCmIDNr9J\ngpHNBksL4CaP9Q4DTuNwZVWrNXmJSOhSwS0iEaWjz5iz3OEwx6ekwZ3XgD3amrwkuNhscEMeLJ5s\njg84jWk2KrpF5HKo4BaRiNHZZxyQ9LyzPTMD/m6msYVQxN31ubCkwBwbdKnoFpHLo5cXEYkInQ7j\nznabj2J75QwV2zK6kklQWmiODbrg/8q9R0mKiPiilxgRCXtdDuPOdmuvOT7962I7ymZNXhI65uXA\nskJzbMBpbKk812FJSiISQlRwi0hY6x2AzUe8i+1p6XD7dBXbcumu81F09zuNkYFajiMiF6KCW0TC\nVv+gcQey2WP035Q0uENtJHIFrsvxPkjZN2j8Une+x/dzRET0ciMiYWnQCVvK4ZzHBsmCVB2QlKtz\nfS7ckGuO9QzA/xyGtl7fzxGRyKaXHBEJO04X/PUkVHksKZk0Hr6hYlvGwMJ87+U4Xf1Ge0lPvzU5\niUjw0suOiIQVlws+Pg0nm83xCUnGnO0YzdmWMTC0HGd2ljne2msU3Y5Ba/ISkeCkgltEwsqeajjc\nYI6lxhsbJOO0QVLGkM0Gt041Rku6a+iC98uNtiYREVDBLSJh5FA97Dtrjo2LhbuLINFuTU4S3mw2\nWDEd8lPM8ao22H7KeMdFREQFt4iEhYoWo5XEXXwM3FUE4+OsyUkiQ3QUfGMWZCWZ4+VN8NkZa3IS\nkeCigltEQt65Dth6HNxvJsZEwT9cA+kJlqUlESQ2Gv7hWkiJN8cP1MGXddbkJCLBQwW3iIS0tl74\n8zFj698QG7BqBkwcb1laEoES7XDXtd7tSzsr4fR5S1ISkSChgltEQlZPP2w+asxAdnfzFJiSbk1O\nEtlSvj6ga/eYhvOXE9DQ6fs5IhL+VHCLSEgadML/lXsvGrk+F4onWpOTCBgjKP9upvFOy5ABp/FO\nTEefZWmJiIVUcItIyHG54MPTUNdhjs/KhEX51uQk4q4gFW6Zao519cN7x8Ax4Ps5IhK+VHCLSMjZ\nXwvHGs2x3GRYPs0Y0yYSDOZkw3yPbZTN3fD+CWMbqohEDhXcIhJSTp2HXVXmWEq8VrZLcFpSANM8\nzhNUtWpcoEik0cuTiISMxi746wlzLC7aGP8Xr8U2EoRsNlg53XtG95d1cKTB93NEJPyo4BaRkNDl\nGGX830xI06xtCWL2aLjzGmPrqbuPTkNtuzU5iUhgqeAWkaA36IT3j0Onwxy/eQpMTrUmJ5HLkRQL\nfz/LWMg0xOmCLcc1uUQkEqjgFpGgt6PSeyLJdRM1/k9CS9Y4WDHNHOvpN9656R+0JicRCQwV3CIS\n1A7VGx/u8lNgaaEl6YhclRmZcEOuOdbUDR+cMsZdikh4UsEtIkGrrgM+qTDHkuOMte1RGv8nIWph\nPkxNM8dONhvjLkUkPKngFpGg1OmALeXmecUxUUYfrCaSSCiz2WDlDMhINMd3V0F1mzU5iYh/qeAW\nkaAz6DSK7e5+c3zldMhM8v0ckVASG2388hgXPRJzAVt1iFIkLKngFpGgs7MS6jvNsQW5MD3DknRE\n/CIlHu6YYY71Dhi/bLqPvxSR0KeCW0SCSnkjlHkckixINfpeRcJNQRrc5PHvdkOX99kFEQltKrhF\nJGic7zaWgbhLjoPbdUhSwtiCXJjicYjySIP3dB4RCV0BLbjLy8u58847ue2220a95t1332XWrFkU\nFxebPvbv3w+Ay+XipZdeYsWKFSxYsIAHH3yQEydOjPr1RCQ09A8ay2363d5Kj7bBN2ZBfIx1eYn4\n29D699R4c/yTCmjssiYnERlbASu4t2zZwiOPPEJBQcFFr83NzaWsrMz0UVJSAsCGDRt49913efnl\nl9mxYwclJSWsXbuWvj6dMhEJVS6XcWf7fI85vmwKTNAhSYkAcTHGL5d2j02U75eDY8C6vERkbASs\n4O7u7ubtt99m0aJFV/V1Nm7cyEMPPcSsWbNITEzkscceo6Ojg507d45RpiISaIcboLzJHJuVCbOz\nrMlHxAoZiXCbxybKtj7YflpLcURCXcAK7nvvvZdJkyZd0rVdXV18//vfZ+HChdx666288847APT2\n9nLy5EmKioqGr7Xb7cycOZOysjK/5C0i/tXQ6X1ALD0Bbp1qvNUuEklmZsKcbHPsZLP6uUVCXdB1\nRqanpzNr1iy+973vMWfOHD766COeeuopsrOzueaaa3C5XKSkpJiek5KSQktLi0UZi8iV6h0w+rbd\nl9vYo+DvZoE9evTniYSz0kI412GsfB+yoxImjleLlUioCropJbfccgtvvvkmJSUlxMbGcscdd7By\n5Uo2b948fI1L762JhDyXCz44Ce0exy9um2bc4RaJVDFR8HczffRzH1c/t0ioCrqC25fc3FwaGhpI\nTU0lKiqK1tZW0+dbW1tJT0+3KDsRuRIH6qDC442p4mzjLXWRSJea4KOfuxc+VD+3SEgKuoJ748aN\nbNmyxRQ7deoU+fn5xMXFMWPGDFO/tsPh4NixY8ybNy/QqYrIFarvhN1V5lhWkvFWuogYfPVzn1A/\nt0hICoqCe9WqVezduxcwCuif//znlJWV0d/fz5///Gd27NjBAw88AMCaNWt46623OH78ON3d3axf\nv56srCyWLFli5bcgIpfIMQh/OWHu246LNt5Cjw6Kn0giwaO0ADITzbGdlZrPLRJqAnZo8o477qC2\nthan08nAwADFxcUAbN26lYqKCrq7jdMhDz74IF1dXTzxxBM0NjaSl5fHyy+/zNy5cwFYvXo1TU1N\nPPzww7S3t1NSUsKrr76K3W4P1LciIlfhkwrjrXF3y6dBcrzv60UiWUw0rJoJbx8cWQo1+HU/97eK\nITboRh+IiC82VxifQKypqWH58uVs376dvLw8q9MRiXjHm4y72+7mZBsjAEVkdL7+27lmgrGhUkSs\nd7GaU2/gikhAtPUa2yTdpSfA0osvnxWJeL76uY81GoW4iAQ/Fdwi4neDTvjrCaN/e0i0De6YoXnb\nIpeqtMB7ZOZHp6Gjz/f1IhI8VHCLiN99XgPnOs2xJQWQqSUeIpcsJtr4JTXKbQOrY9D4ZdYZts2h\nIuFBBbeI+FVNG+w7a45NSYO5E63JRySUZSYZv6y6q+2A/bXW5CMil0YFt4j4TU8//PWkOZZkN6aS\n2Gy+nyMiF3bdRJicYo7trTbm24tIcFLBLSJ+4XLBh6egy2GOr5wBCZriKXLFbDZYMR3i3UYCOl3e\n5yREJHio4BYRvzjaCKc9VrdfPwnyU3xfLyKXLinWeKfIXWuvsRRHRIKPCm4RGXNtvbCjwhzLHgcL\n863JRyQcTU33HhV4pAFONluTj4iMTgW3iIwppws+ODmyFQ8gJgpun67V7SJjbWkBpHpsaf3oNHRq\nVKBIUNHLn4iMqS/rjKkJ7pYWQGqC7+tF5MrZfYwK7B2AD04Z5yhEJDio4BaRMdPcDburzLHJKd5v\ne4vI2MkaBzd5tGtVt8GhemvyERFvKrhFZEwMbZN0X8ARFw3Lp2sEoIi/lUyC3GRz7NMzxnkKEbGe\nCm4RGROf10BTtzl2y1QYF2tNPiKRxGaDFdPA7vaqPuCEbSe1hVIkGKjgFpGrVtcBf/PYJjkjA2Zm\nWpOPSCRKjofSQnOsrgO+qrMkHRFxo4JbRK5K/6BxF839JlqSHW6ZYllKIhGrKAsKUs2x3VVwvtv3\n9SISGCq4ReSq+OoTXT4d4rVNUiTgbDa4bZpxfmLIoEutJSJWU8EtIlesqtV7EkJxtvcdNhEJnHGx\ncLPHO0wNXd5tXyISOCq4ReSKOAZg+ylzLDUelhRYk4+IjJiZCdPSzbHPa6Cxy5p8RCKdCm4RuSKf\nnoFOx8hjG7BiurGIQ0SsZbPBrVMhIWYk5vy6tWTQOfrzRMQ/VHCLyGWrboPDDebYvBzIGW9NPiLi\nLcEOt04zx5q7YW+1NfmIRDIV3CJyWRyD8KGPVhLPTXciYr1p6TDLYzzn/lqo77QmH5FIpYJbRC7L\nrjPQ3meOrZgOMWolEQlKy6ZAktsCKhewXa0lIgGlgltELllNG5R5TCVRK4lIcIuPgeVTzbHmHtin\nqSUiAaOCW0QuSf+g91SSlDi1koiEgoI0uGaCObbvLDRpaolIQKjgFpFLsrvKu5VkuaaSiISM0gJI\ndFtI5XQZv0RrIY6I/6ngFpGLqm2Hr86ZY3MnQm6yNfmIyOWLt/teiPNlrTX5iEQSFdwickG+WkmS\n42DxZGvyEZErNz3DeyHOnmpo6bEmH5FIoYJbRC5obzW09ppjy6eplUQkVN08xThIOWTw69YSl1pL\nRPxGBbeIjKqhE76sM8eKsyEvxZp8ROTqJcVCaaE5VtfhPYFIRMaOCm4R8WnQ+fVdL7fY+FhYXGBZ\nSiIyRmZlQkGqObbrDLT3+r5eRK6OCm4R8enLOmjqNsdumQqxaiURCXk2G9w61dwa1u+ED0+rtUTE\nH1Rwi4iXtl6jd9vdzAwoTLMmHxEZe+PjYInH4efqNjjWaE0+IuFMBbeImLhcxl2uQbe7XPExUDpl\n9OeISGiak+093vPTM9DTb00+IuFKBbeImBxrNFa4u1vqsTBDRMKDzQa3TYVo20isdwB2VlqWkkhY\nUsEtIsO6+427W+7yU7xXQotI+EhNgBvzzbHyJqhqtSYfkXAU0IK7vLycO++8k9tuu+2C123bto27\n776b+fPns3LlSl5//fXhzz3zzDNce+21FBcXD3/Mnz/f36mLRISdFcbdrSExUcbBKptt9OeISOib\nnwMZiebYR6eNxVcicvUCVnBv2bKFRx55hIKCC88UO3jwIE899RSPPvooX3zxBevWreM3v/kNW7du\nHb7mrrvuoqysbPjjwIED/k5fJOxVtsDxZnNsYR6kxFuTj4gETnSU0Vrirr0PvqixJh+RcBOwgru7\nu5u3336bRYsWXfC61tZW1q5dy6pVq4iJiWHBggVcf/317Nu3L0CZikQexyB8fNocm5AE8yZZk4+I\nBN7E8TB3ojl2oA6auqzJRyScBKzgvvfee5k06eKv3suWLePxxx8ffuxyuaivrycrK2s4Vl5ezurV\nqykpKeHuu+9m//79fslZJFLsqYIOx8hjG8bdrii1kohElEX5xibKIc6vpxY5NZtb5KoE/aHJ1157\njdbWVu6//34A8vPzKSws5MUXX2Tnzp2UlpbyyCOP0NzcfJGvJCK+1HfCwXPm2LwcyBpnTT4iYp3Y\nGLjZYwRofScc0tp3kasS1AX3yy+/zO9+9ztee+01UlONHbSPPfYY69evJy8vj6SkJJ588kkSExPZ\ntm2bxdmKhB6ny2glcb95lRwHC/NHfYqIhLlp6TDVY8nVriro7LMmH5FwEJQFt8vl4qc//Sl/+tOf\n2LBhA0VFRaNeGx0dTU5ODg0NDQHMUCQ8lJ2DBo/+zFummNc9i0jkudnj50D/IHxSaVk6IiEvKAvu\nX/7yl3z55Zf84Q9/YNq0acNxl8vFunXrOHbs2HDM4XBQVVVFfr5uyYlcji4H7PFY3z49Awq0vl0k\n4o2LM/q53Z0+D6fOW5OPSKgLioJ71apV7N27F4D9+/ezadMmfvvb35KZmWm6zmazUVNTw7PPPkt9\nfT1dXV288MIL2O12br/9ditSFwlZOyuN6SRD7NFQWmhVNiISbIonQrbHWY4dFeafGyJyaWIC9Rfd\ncccd1NbW4nQ6GRgYoLi4GICtW7dSUVFBd3c3AJs2baK7u5uVK1eann/DDTfwu9/9jueee45f/epX\n3HPPPXR2djJ37lzeeOMNkpKSAvWtiIS8qlY44XHOeFE+jIv1fb2IRJ6or9e+v102MqWk02HM5l5y\n4ZUaIuLB5nK5wnbYT01NDcuXL2f79u3k5eVZnY5IUBgYhA1fQZvbAagJSXB/scYAioi3z87A/tqR\nx1E2+NZc782UIpHsYjVnULSUiEjg7Ks1F9tgrG9XsS0ivtyQZ373y+ky1r6H7+06kbGnglskgrT0\nwN/OmmPF2d59miIiQ2KjYVmhOVbXAccaLUlHJCSp4BaJEC4XfFxh3hiXaIdFk63LSURCw9R0KEw1\nxz47A7391uQjEmpUcItEiONNUNNmjpUWQlzAjk6LSKiy2WDZFIh2az3rGYDd1aM/R0RGqOAWiQB9\nA7DzjDmWnwIzMqzJR0RCT0q80c/t7lA9nOuwJh+RUKKCWyQC7K6CHre3fqNsxiY5mw5KishlKJkE\nqfHmmGermoh4U8EtEubOdUBZvTm2IBfSEqzJR0RCV3SU8cu6u8YuKDtnTT4ioUIFt0gYc359UNJd\nSjxcn2tNPiIS+ianerej7amGLoc1+YiEAhXcImHs4Dnj7pO7W6ZAjP7LF5GrUFoI9uiRx45B+PTM\nqJeLRDy97IqEqS6HcdfJ3YwM4+6UiMjVSIqFRfnm2PEmqG7zfb1IpFPBLRKmPjsD/YMjj2OjjbtS\nIiJjoXgiTEgyxz4+DYNOa/IRCWYquEXCUG07lDeZYwvzjbtSIiJjIcpmtKi5a+2F/bXW5CMSzFRw\ni4QZXwclMxJh7kRr8hGR8DVxPMzOMse+qIG2XmvyEQlWKrhFwkzZOWjuNsdunmLcjRIRGWuLJ0O8\n28baQRfsrLQsHZGgpIJbJIx093sflJyVCbnJ1uQjIuEv3g5LC8yxihaobLEmH5FgpIJbJIzsOmOM\n5xpij4LFBaNfLyIyFq6ZADnjzbGdlTpAKTJEBbdImKjrgKON5tiN+TBOByVFxM9sNu8NlK29cKDO\nmnxEgo0KbpEw4HTBJx4HJdMS4DodlBSRAJmQBHOyzbF9NdDZZ00+IsFEBbdIGDhc771R8uYpEK3/\nwkUkgG7KNx+g7HdqA6UIqOAWCXk9/bC7yhybngH5KdbkIyKRK8EOiyabYyeaoUYbKCXCqeAWCXG7\nq6DP46BkqQ5KiohFirK8N1DuqNABSolsKrhFQlh9JxxuMMduyINxcdbkIyIS5eMAZXMPlNVbk49I\nMFDBLRKiXD4OSqbGw7wca/IRERmSM94YFehubzV0O6zJR8RqKrhFQtSRBuMOtzsdlBSRYLF4MsRG\njzx2DMKuqtGvFwlnemkWCUG9/d4vXNPSYXKqNfmIiHhKioUb88yxo41wrsOafESspIJbJATtqYbe\ngZHHMVHiR0ZvAAAgAElEQVSwtNCydEREfJo7EdITzLFPKozdASKRRAW3SIhp7oZDHoePFuRCsg5K\nikiQiY6CZR4HKBu6jJY4kUiiglskhLhcsLMS3G8OJcfB/ElWZSQicmH5KcZuAHe7q4zWOJFIoYJb\nJIRUtEC1xwKJpQVGS4mISLDy/DnVO2C0xolECr1Mi4SIQSd8WmmO5SXD1HRL0hERuWTj44zWN3eH\n6qGxy5p8RAJNBbdIiPiqDtr6Rh7bgNJCsNmsykhE5NKVTIIUt7MmLowNlC4doJQIoIJbJAR0O+Dz\ns+bYnGzITPJ9vYhIsImOglKPA5S1HXDyvDX5iASSCm6RELC7GvoHRx7HRcPCfOvyERG5EoWp3vsC\nPqs0/3wTCUcquEWCXEOn9witG/MhwW5NPiIiV8pmM1rhotxa4ToccKDWspREAkIFt0gQc7lgR6U5\nlpYAxdmWpCMictXSE4yFOO7+Vgsdfb6vFwkHAS24y8vLufPOO7ntttsueN3WrVu56667mD9/Pv/4\nj//IX//61+HPuVwuXnrpJVasWMGCBQt48MEHOXHihL9TF7HEyWao81iDvLTA6IUUEQlVN+ZBQszI\n4wEnfHbGunxE/C1gL9tbtmzhkUceoaCg4ILXHTt2jKeffpof/OAH7NmzhyeeeIIf/ehHHD9+HIAN\nGzbw7rvv8vLLL7Njxw5KSkpYu3YtfX361VjCy8Cg9wtQQSoUplmTj4jIWImLgZsmm2MnmuFsuzX5\niPhbwAru7u5u3n77bRYtWnTB69555x2WLFnCihUriIuLY/ny5SxatIg//vGPAGzcuJGHHnqIWbNm\nkZiYyGOPPUZHRwc7d+4MxLchEjD764zexiFRX/c+ioiEg6IsmOAxaWlnJTg1JlDCUMAK7nvvvZdJ\nky6+f/rw4cPMnj3bFCsqKqKsrIze3l5OnjxJUVHR8OfsdjszZ86krKxszHMWsUpnH/zNYwzg3IlG\n/7aISDjwdROhsQuONvi8XCSkBV0naGtrK8nJyaZYSkoKLS0ttLW14XK5SElJ8fl5kXCxq8roaRwS\nH2P0PIqIhJPcZJiRYY7troK+AWvyEfGXoCu4wTgYeTWfFwlldR1Q3mSOLZps9DyKiISbJQUQ41aN\n9AzA5zXW5SPiD0FXcKelpdHa2mqKtba2kpGRQWpqKlFRUT4/n56eHsg0RfzC5TJWHbvLSDR6HUVE\nwtH4OGPtu7uD56Clx5p8RPzhsgpuh8PByZMn+fzzz9m7dy8nTpzA4XBc/ImXYc6cORw6dMgUKysr\n47rrriMuLo4ZM2aY+rUdDgfHjh1j3rx5Y5qHiBWONUFDlzm2rNC8JEJEJNyUTIJxsSOPnS7jAKVI\nuLikgvuDDz7gO9/5DjfccAN33nknDz74IA899BD/8A//wA033MB3vvMdPvjggytOYtWqVezduxeA\nb33rW+zdu5dt27bhcDh4//332bdvH9/61rcAWLNmDW+99RbHjx+nu7ub9evXk5WVxZIlS6747xcJ\nBo5B2O0xBnBaOuSl+L5eRCRc2KON1hJ3Z1qhUsezJExcsCu0urqaH/7wh5w9e5a77rqL7373u8yY\nMYO0tDRsNhvnz5/nxIkTfP755/zsZz/jt7/9LS+++CJ5ed6nu+644w5qa2txOp0MDAxQXFwMGEtu\nKioq6O7uBmD69OmsX7+eF198kX/5l3+hsLCQ//zP/xye37169Wqampp4+OGHaW9vp6SkhFdffRW7\nXXuuJbT97Sx09Y88jrJ5vwCJiISrGRlQdg5q3ZZ97ayE/BQt+5LQZ3Nd4ATi0qVL+d73vscDDzxA\nbGzsaJcBRmvHhg0beP311/n000/HPNErUVNTw/Lly9m+fbvPXwJEgkVbL/z3lzDo9l/j9bmwePLo\nzxERCTeNXfCHg+bY0gKYf/GpwiKWuljNecE73G+++SZTp069pL8oNjaW73znOyxbtuzKMhWJYJ+d\nMRfbiXZYkGtdPiIiVpiQBLOz4LDbLO7Pa2BWJiRe+L6fSFC74Js07sX2Qw89xNatW31ed9111/l8\njohcXE0bnDpvji2eDLHR1uQjImKlRR4//xyDsLvaunxExsIld0V98cUXPPvss/zsZz/zmkyiudgi\nV8bXSfysJLhmgiXpiIhYLsHuvejrSAM0dFqTj8hYuOSC22638z//8z+cPHmS++67j8rKyuHP2Wya\nWSZyJY40QFO3ObZsCug/KRGJZHMnQlqCObaj0thVIBKKLuvcb3Z2Nm+++Sa33XYb9957L5s3bwZ0\nh1vkSvQNwJ4qc2xmJuSMtyYfEZFgER0FpR5Tmuo64ESzNfmIXK3LXhYdFRXFE088wY033si//uu/\nsmfPHn/kJRL2Pq8xVhgPiYnSVBIRkSEFaVCYCpVuy6U/OwNT0oy53SKh5JLvcOfk5JgeL1q0iM2b\nN9PY2Eh/f/8ozxIRX1p6jNXF7q6fZKw4FhERw9JC86bdTgfsr7UsHZErdsGCu76+fvh/+5pQkp6e\nzuuvv8727duHYw0NDV7XiYjZp5XGgckh42I1Z1ZExFNaAlw30RzbXwsdfdbkI3KlLlhwf/Ob3+S9\n99676BeZNMmoFN577z2++c1vjk1mImHqTKv5LVIwNkrqLVIREW835EGCWwPsgBN2nbEuH5ErccEe\n7ldffZWnnnqK3/zmN9x///3cdNNNzJw5c3iNen9/PydOnGDPnj388Y9/xOl08uqrrwYkcZFQNOj0\nHgOYM95YaSwiIt7iYuCmyfDR6ZHY8WaY26FD5hI6Llhwz549m/fee49NmzaxceNGnn/+eWw2G/Hx\n8dhsNnp6enC5XMyYMYMHH3yQb37zmxddAS8SyQ7VG/3b7pYVagygiMiFFGVB2TnzGNUdFXB/sX5+\nSmi46JSS2NhYvv3tb/Ptb3+bpqYmDh48yKlTp7DZbMyYMYOioiImTNCWDpGL6emHvR7b0oomQNY4\na/IREQkVUTbj5sS7R0ZiDV1wrBGuzbIsLZFLdsljARsaGvi3f/s3du/ezcCAMcssJiaG0tJSfv7z\nn5OZmem3JEXCwd5q6BsceWyPNt4mFRGRi8tNgenpcPL8SGxXFUzLMK+CFwlGlzwW8Mknnxzu0X7/\n/ffZsmULr7zyCg6HgyeeeMKfOYqEvKYuo53E3Q25kKQOLBGRS7akAKLdWki6+2HfWevyEblUl3yH\n+9ChQ+zatYtx40be/546dSpz585l2bJlfklOJBy4XLDzDLjvY02Jg3k5oz5FRER8SI43Rqi6F9kH\namF2FqTEW5eXyMVc8h3ugoICurq6vOIOh4PJk/W+uMhoKlqgps0cW1porC4WEZHLc30uJNlHHjtd\nxgZKkWB2yXe4f/CDH/CjH/2IBx54gClTpjA4OEhVVRVvv/023/3udzl58uTwtdOnT/dLsiKhZtBp\nLLlxl5dirCYWEZHLFxsNiwrgg5Gyg1PnjRsbeSnW5SVyIZdccP/zP/8zAF988YXX5/bu3YvNZsPl\ncmGz2Th69OjYZSgSwr6sgza3jWg2oLRAY6xERK7GNZnGmMD6zpHYjkr41lzzKniRYHHJBbf7+nYR\nubguB3zhcZhnTjZkJlmTj4hIuLDZoLQQNh0aiTV3w5EG4+esSLC55II7NzfXn3mIhJ09VdDvNgYw\nLhoW5luXj4hIOMkZD7MyobxpJLa7ytjcG3fJ1Y1IYOjYlogfNHTCkUZzbGE+JNh9Xy8iIpdv8WSI\ncatkegfg8xrr8hEZjQpukTHmchm9hO7SEvQ2p4jIWBsXB9dPMscOnoOWHmvyERmNCm6RMXaiGeo6\nzLHSQo0BFBHxh/mTYLzbEjGny3s6lIjVVAKIjKH+QdjlMQ+2MBUKUq3JR0Qk3NmjYXGBOVbZCmda\nrMlHxBcV3CJj6EAtdDhGHkfZjCU3IiLiPzMyjEOU7naeMXYhiAQDFdwiY6SzD/5Wa47NnWj0b4uI\niP/YbLCs0Bxr6YFD9ZakI+JFBbfIGNlVBQNud1PiY+DGPOvyERGJJFnjoGiCOba3Gnr6rclHxJ0K\nbpExUNdhngULsGiyZsGKiATSTZONnu4hfYNG0S1iNRXcIlfJ5YIdFeZYZiIUZVmTj4hIpEqKhRs8\n9vQdqje2UIpYSQW3yFU61ggNXeZYaaFxYFJERAJrXg6kxI08dmHsRnC5rMpIRAW3yFVxDBqrhN1N\nS4e8FGvyERGJdNFRsKTQHKtpgwqNCRQLqeAWuQr7zkKX24GcaBssKRj9ehER8b+pad43Pj6t1JhA\nsY4KbpEr1NYLX3qMAZw/CVLirclHREQMNhuUFoB7Z19bH3xVZ1lKEuFUcItcoc/OwKBbT2CSHa7P\nHf16EREJnMwkmJNtjn1+Frodvq8X8ScV3CJXoKYNTp03xxYVQGy07+tFRCTwFuZDnNvP5f5B2K0x\ngWIBFdwil8npMk68u8seB9dkWpKOiIiMIsEON+abY0caoLHL9/Ui/hLQtRx1dXU8++yzHDhwgPj4\neJYvX84zzzxDbGys6bqf/OQnbN682RRzOp2UlJTw1ltv8cwzz7B582ZiYkbSj4mJ4cCBAwH5PiSy\nHWnwnulaWmj0DIqISHApzjZmcbf0jMR2VMA9s/VzWwInoHe4H3/8cdLS0ti2bRsbNmzgwIEDvPTS\nS17X/eIXv6CsrGz44+DBgxQXF3PnnXcOX3PXXXeZrlGxLYHQN+A9BnBWJuSMtyYfERG5sOgoWOox\nPaq2A042W5OPRKaAFdxlZWUcOXKEp59+muTkZHJzc1m7di3vvPMOTueF5/Rs2rSJ/v5+7rvvvgBl\nK+Lb5zXQOzDyOCYKFk+2Lh8REbm4wjQoSDXHPjsDA4PW5CORJ2AF9+HDh8nJySE9PX04Nnv2bNra\n2qiqqhr1eT09Paxfv54f//jHREWNpFteXs7q1aspKSnh7rvvZv/+/X7NX6SlBw6eM8euz4Vxcb6v\nFxGR4OG5AbjDAQc0JlACJGAFd2trK8nJyaZYSooxlb6lZfT1Txs3bmTatGksWLBgOJafn09hYSEv\nvvgiO3fupLS0lEceeYTmZr0/JP7zaaVxYHLI+FgoybEsHRERuQxpCTB3ojm27yx09lmTj0SWgPZw\nu1yui1/kxul08vvf/56HH37YFH/sscdYv349eXl5JCUl8eSTT5KYmMi2bdvGMl2RYWdaoLLVHFtS\nADEaAygiEjJuyIN4t3ERA07YNfqb7CJjJmAFd3p6Oq2t5opl6HFGRobP5+zbt4/Ozk6WLFlywa8d\nHR1NTk4ODQ0NY5OsiJtBJ+w8Y45NGg/Tff9rKyIiQSo+Bm7yGBNY3gTnOqzJRyJHwAruOXPmUF9f\nT2Nj43Ds4MGDZGRkkJ+f7/M5H3zwATfddJNpbKDL5WLdunUcO3ZsOOZwOKiqqhr164hcjTKPcVKg\nMYAiIqFqdjZkJJpjOyrhMt+EF7ksASu4i4qKmDdvHs8//zwdHR1UV1fzyiuvsGbNGmw2G6tWrWLv\n3r2m5xw5coS8vDxTzGazUVNTw7PPPkt9fT1dXV288MIL2O12br/99kB9OxIhevphr8dWsqIsyBpn\nTT4iInJ1omzGTRN39Z3GnW4RfwloD/evf/1r2tvbKS0t5b777mPZsmU8+uijAFRUVNDdbd4m0tjY\nyIQJE7y+znPPPUdhYSH33HMPixcv5ujRo7zxxhskJSUF5PuQyLGnGhxuY6Nio2GR3kgREQlp+Skw\nNd0c23XG/PNeZCwFdNNkdnY2//Vf/+Xzc+Xl5V6xv/zlLz6vTU1NZd26dWOam4inxi5jO5m7G/Ig\nMdb39SIiEjqWFkBly8j0qa5++NtZWKTdCuIHAb3DLRIqXC5j9a+71Hi4bqLv60VEJLSkxMN8j9Gu\nB2qhvdeafCS8qeAW8eFks7H6111pobEiWEREwsOCPEi0jzwedMFnGhMofqDyQcRD/yB86jEGsCDV\nWA0sIiLhIzbau4XkZDOcbbcmHwlfKrhFPOyvhU7HyGNfJ9pFRCQ8XDsBsjxmLuyoMG8WFrlaKrhF\n3HT0GYdm3F030VgJLCIi4cfm46ZKUzcc1S49GUMquEXcfHrG6OEbkmA3JpOIiEj4mpQMMzy2B++u\nhr4Ba/KR8KOCW+RrZ9uN3j13i/IhLqDDM0VExApLCiDGrSrq6YcvaqzLR8KLCm4RjF49zzGAWUnG\nVkkREQl/4+OgZJI59tU5aO2xJh8JLyq4RYAjDUbPnrtlU4zePhERiQwlk2Cc23Izp8t7apXIlVDB\nLRGvdwB2e8xdnZkJOeOtyUdERKxhj4bFHmMCK1qgqtWafCR8qOCWiPd5tVF0D4mJgiVa7SsiEpF8\n3XDZWakxgXJ1VHBLRDvfDQfPmWMLcmFcnDX5iIiItXyNCTzfA4fqLUlHwoQKbolYLhfsqAT3mxbJ\ncTA/x6qMREQkGGSPg2smmGN7q6G335p8JPSp4JaIVdEC1W3m2JICiIm2Jh8REQkeiyeD3a1K6h2A\nvRoTKFdIBbdEpEEnfFppjuUlw7R0S9IREZEgkxQLCzwWn5WdM1oRRS6XCm6JSF/WQVvfyGMbRs+e\nxgCKiMiQeTlGq+EQF8YBSpcOUMplUsEtEafL4b09bE42ZCZZk4+IiASnmCij1dBdVRtUakygXCYV\n3BJxPjsD/c6Rx3HRsDDfunxERCR4TUuH3GRzbGcFDDh9Xy/iiwpuiSi17VDeZI7dNBkS7NbkIyIi\nwW1oTKB7x2FbHxyotSojCUUquCViOF3wcYU5lpFotJOIiIiMZkKS92vFvrPQ0ef7ehFPKrglYpSd\ng2aP0+U3T4EoHZQUEZGLuCkf4mNGHg84jQOUIpdCBbdEhO5+2FNtjs3K9O7LExER8SXeDosmm2On\nzkOVDlDKJVDBLRFh1xlwDI48tkd7nzwXERG5kNlZkOUx0eqTCmO3g8iFqOCWsHeuA442mmML84yl\nBiIiIpfKZoNbpppjrb3GbgeRC1HBLWHN10HJ9ASYO9GafEREJLRlj4OiLHPsixro1AFKuQAV3BLW\nDtdDY5c5dvMUiNa/+SIicoUWTzZ2OAzpd8KnZ6zLR4Kfyg4JWz0+DkrOyIC8FGvyERGR8JDg4wDl\niWaoabMmHwl+KrglbO2ugt6Bkcf2KFiqg5IiIjIGZmcb87nd6QCljEYFt4Sl+k443GCO3ZAH4+Ks\nyUdERMJLlM1oUXR3vgcOnrMmHwluKrgl7Lhcxl0Gd6nxMC/HmnxERCQ85YyHayaYY3troMthTT4S\nvFRwS9g50mDc4Xang5IiIuIPSyZDrPsBykH4TAcoxYNKEAkrvf2wq8ocm5YOk1OtyUdERMJbYqyx\n9t1deRPUtluTjwQnFdwSVvZUmw9KxkTB0kLL0hERkQhQPBEyEs2xjyuMXRAioIJbwkhjFxyqN8cW\n5EKyDkqKiIgf+TpA2dwNZTpAKV9TwS1hwemCD0+B+82ElDiYP8mylEREJILkJsOsTHNsdzV06gCl\nEOCCu66ujkcffZSFCxdy88038x//8R84HN7/Jr777rvMmjWL4uJi08f+/fsBcLlcvPTSS6xYsYIF\nCxbw4IMPcuLEiUB+KxJkDtdDg8dGydIpRkuJiIhIICwuALvHAcpPKy1LR4JIQMuRxx9/nLS0NLZt\n28aGDRs4cOAAL730ks9rc3NzKSsrM32UlJQAsGHDBt59911efvllduzYQUlJCWvXrqWvry+Q344E\niW6H74OSU9KsyUdERCLTuFhY5HGA8kQzVLVak48Ej4AV3GVlZRw5coSnn36a5ORkcnNzWbt2Le+8\n8w5O5+WtZdq4cSMPPfQQs2bNIjExkccee4yOjg527tzpp+wlmH16BhyDI4/tUbCs0LJ0REQkghVP\n9N5A+XEFDGgDZUQLWMF9+PBhcnJySE9PH47Nnj2btrY2qqqqvK7v6uri+9//PgsXLuTWW2/lnXfe\nAaC3t5eTJ09SVFQ0fK3dbmfmzJmUlZX5/xuRoFLdZoxfcrcwXxslRUTEGlE2uHWqOdbWC387a00+\nEhxiAvUXtba2kpycbIqlpKQA0NLSQmFh4XA8PT2dWbNm8b3vfY85c+bw0Ucf8dRTT5Gdnc0111yD\ny+Uafq7712ppafH79yHBY9AJH582xzIT4TptlBQREQtlj4PibChzm5y17yzMzIS0BOvyEusEtIfb\n5bq0gZS33HILb775JiUlJcTGxnLHHXewcuVKNm/efNlfS8LX/lpo7TXHbplq3F0QERGx0qLJkGgf\neex0wScVoPIlMgWs4E5PT6e11XxqYOhxRkbGRZ+fm5tLQ0MDqampREVF+fxa7u0qEt5ae+CLGnNs\ndhbkjLcmHxEREXdxMbC0wByrbjMOUUrkCVjBPWfOHOrr62lsbByOHTx4kIyMDPLzzUd6N27cyJYt\nW0yxU6dOkZ+fT1xcHDNmzDD1azscDo4dO8a8efP8+01IUHB9fZdg0O0uQUIMLJ5sXU4iIiKeZmZC\nnrkDlp2V0Dfg83IJYwEruIuKipg3bx7PP/88HR0dVFdX88orr7BmzRpsNhurVq1i7969gFFA//zn\nP6esrIz+/n7+/Oc/s2PHDh544AEA1qxZw1tvvcXx48fp7u5m/fr1ZGVlsWTJkkB9O2Khk+ehqs0c\nW1II8Xafl4uIiFjCZoNbpphbHbv7Ybf3rAgJcwE7NAnw61//mn//93+ntLSU+Ph4/umf/olHH30U\ngIqKCrq7uwF48MEH6erq4oknnqCxsZG8vDxefvll5s6dC8Dq1atpamri4Ycfpr29nZKSEl599VXs\ndlVc4c4xADsrzLHcZLgm0/f1IiIiVkpLgOtzzW2QZfVwbZZxuFIig80VxqcPa2pqWL58Odu3bycv\nL8/qdGQMfFIBB8+NPI6ywQNzIT3RupxEREQuZMAJG76ENrf9fBOS4P5iHfQPFxerObX4WkLGuQ5z\nsQ1QMknFtoiIBLeYKLjZYzZ3Yxd8WWdNPhJ4KrglJAw64cNT5lhyHCzItSYfERGRy1GQCjM8hrLt\nrTaW4kj4U8EtIWF/LTT3mGO3TgV7tDX5iIiIXK5lhRDn9ro14ISPTms2dyRQwS1Br8XHzO1rJsDk\nVGvyERERuRKJsbC00ByrboPyJkvSkQBSwS1BzeUyfvt3n7kd72OZgIiISCi4dgLkJZtjOyuhp9+S\ndCRAVHBLUDvcAGfbzbFlhZCgCZAiIhKCbDajJTLabTpJ74BRdEv4UsEtQavLAZ+dMccmpxqbu0RE\nREJVagLcaF6yTXkTnGm1Jh/xPxXcErR2VIBjcORxTBTcOsW4OyAiIhLK5udApsdY249PQ/+g7+sl\ntKnglqB0+ryxwt3dTfmQHG9NPiIiImMpOgpumwbu95Da+2BPtWUpiR+p4Jag4xiAjz3Wt2clwXU5\n1uQjIiLiD9njvF/bvqqD+k5r8hH/UcEtQWdXldG/PcSGcRdA629FRCTcLMyH8XEjj10Yi94GnZal\nJH6ggluCytk2KKs3x0omwYQka/IRERHxp9houGWKOdbUbSx8k/ChgluCRv8gfOCxvj0lDm7MsyYf\nERGRQChM857A9XkNNHVZk4+MPRXcEjR2VRkHRtwtnwYxWt8uIiJhblkhJMSMPHa6jJtQai0JDyq4\nJSicbYOD58yxuRMhN8WafERERAIpwQ63TDXHGrvUWhIuVHCL5Xy1kiTHweLJ1uQjIiJihekZMCPD\nHFNrSXhQwS2W89VKsmIa2NVKIiIiEebmKWotCUcquMVSaiUREREZodaS8KSCWyyjVhIRERFvai0J\nPyq4xTJqJREREfFNrSXhRQW3WEKtJCIiIqNTa0l4UcEtAadWEhERkYtTa0n4UMEtAffpGbWSiIiI\nXApfrSV/PanWklCjglsCqqIFDtWbY2olERER8c1Xa0lzN+yptiYfuTIquCVgevrhQ49WktR4tZKI\niIhcyPQMmJVpju2vhbPt1uQjl08FtwSEy2UU2939IzEbsHK6WklEREQu5uYpMC7WHNt2EvoGrMlH\nLo8KbgmIo41wusUcuyEPJo63Jh8REZFQEhcDK6abYx19sKPSknTkMqngFr9r64UdFeZY9jij4BYR\nEZFLk58C83LMsWONcLLZmnzk0qngFr9yuoy3vPrdTlPHRBmtJFE26/ISEREJRYsmQ3qCOfbRaeh0\nWJOPXBoV3OJX+2qgrsMcW1oAaQm+rxcREZHRxUTB7TPMN616B4ybWy6XdXnJhangFr+pbTcG9Lsr\nSIU52dbkIyIiEg4mJMFN+eZYTZu2UAYzFdziF70D8JcT4P7LdoIdlk8Dm1pJRERErsr8SZCbbI7t\nqYZzHb6vF2up4JYxNzQC0LOfbOU0SIr1/RwRERG5dFE2uH06xHtsofzLCY0KDEYquGXMHW6AU+fN\nsfk5UJBmTT4iIiLhaFyc8c6xu/Y++Pi0+rmDjQpuGVPN3d4jALOSjFPVIiIiMrampsPciebY8WZj\nXKAEj4AW3HV1dTz66KMsXLiQm2++mf/4j//A4fA9x2bbtm3cfffdzJ8/n5UrV/L6668Pf+6ZZ57h\n2muvpbi4ePhj/vz5gfo2ZBT9g/CX4zDo9lu1PQrumAHR+tVORETEL5YUQEaiOfZJBZzvsSYf8RbQ\nMujxxx8nLS2Nbdu2sWHDBg4cOMBLL73kdd3Bgwd56qmnePTRR/niiy9Yt24dv/nNb9i6devwNXfd\ndRdlZWXDHwcOHAjktyI+fFIBzR7/cd88FVI1AlBERMRvYqJg1QzjzyH9TthabtwME+sFrOAuKyvj\nyJEjPP300yQnJ5Obm8vatWt55513cDqdpmtbW1tZu3Ytq1atIiYmhgULFnD99dezb9++QKUrl+lI\ng7G+3d2sTLgm05p8REREIkl6IpQWmmPNPfBxhfq5g0HACu7Dhw+Tk5NDenr6cGz27Nm0tbVRVVVl\nunbZsmU8/vjjw49dLhf19fVkZWUNx8rLy1m9ejUlJSXcfffd7N+/3//fhPjU1GUc0HCXGg+3TNUI\nQJUhFU8AABZBSURBVBERkUCZnQUzPW50HWs0boqJtQJWcLe2tpKcbB4YmZKSAkBLS8sFn/vaa6/R\n2trK/fffD0B+fj6FhYW8+OKL7Ny5k9LSUh555BGam5v9k7yMyjEA73v0bcdEwTdmQWy0dXmJiIhE\nGpsNbp3qvc35kwpo7LImJzEEtIfbdQXvabz88sv87ne/47XXXiM1NRWAxx57jPXr15OXl0dSUhJP\nPvkkiYmJbNu2baxTlgtwuWD7aWjtNcdvmep9eENERET8LzYavjHT3M896DJujmk+t3UCVnCnp6fT\n2tpqig09zsjI8Lre5XLx05/+lD/96U9s2LCBoqKiUb92dHQ0OTk5NDToPZNAOngOTnq8qVCUBddO\nsCYfERERMfq5b51qjrX1wvZT6ue2SsAK7jlz5lBfX09j48jJuoMHD5KRkUF+fr7X9b/85S/58ssv\n+cMf/sC0aSNT3V0uF+vWrePYsWPDMYfDQVVVlc+vI/5xth0+PWOOZSbCzYWWpCMiIiJurpkAc7LN\nsVPnYX+tNflEuoAV3EVFRcybN4/nn3+ejo4OqqureeWVV1izZg02m41Vq1axd+9eAPbv38+mTZv4\n7W9/S2amufvfZrNRU1PDs88+S319PV1dXbzwwgvY7XZuv/32QH07Ea2jD94vN1bIDomNNvq2Y9S3\nLSIiEhRKC2FCkjm2qwrOXPjonPhBQHu4f/3rX9Pe3k5paSn33Xcfy5Yt49FHHwWgoqKC7u5uADZt\n2kR3dzcrV640Lbd5+OGHAXjuuecoLCzknnvuYfHixRw9epQ33niDpKSkUf9uGRsDg7ClHHo8+sBW\nTIeUeGtyEhEREW8xUfB3MyHO42bYX05Aq5biBJTNdSUnGUNETU0Ny5cvZ/v27eTl5VmdTshzuWDb\nSShvMsdvzIOF6uYREREJSlWt8L9Hwb3gS0+A+4o1UWysXKzm1MJtuWRfnfMutqekGQW3iIiIBKfJ\nqbC4wBw73wPbTugQZaCo4JZLUt0Gn1aaY2kJcPt0LbcREREJdvNzvJfinG6BL2qsySfSqOCWizrf\nYxySdP8lODYa/n4WxMZYlpaIiIhcIpsNbpvqfYhybw0cb/L9HBk7Krjlgnr64b2j0Ddojt8+w3uT\nlYiIiAQv+9cTxeI9bpZ9cBLqOqzJKVKo4JZRDTjh/8r/f3v3Hhx1ee9x/JN7FmFzI4EQIlEg6CaU\nBLE0UhTCLefUC6FUM4NnOrWdgdr0gpWpY4di2xlrBx1Go4OCtZ2p3GLQhqNAyVHkQOcQSUGzAVdE\nsCFcwoLZJLDEkOyeP35NZJMg2Q17yeb9mtnRPL8n8bs+G/jsb5+L1PKlZ3vBzcbcbQAAMLiY44yT\nKCOvmg7a6ZbesRmH48A/CNzok9vd9zteS6p0x5jg1AQAAAYuI8GYXnK1yx3Sf9ukNo5/9wsCN/pU\nfVL6tMex7WPN0qxbWSQJAMBgd3uadGeGZ1vTv9dsdbqCU1M4I3CjlyPnpAOnPNuSTNJ/TJKieMUA\nABAWpmdKE1M82xpapN3H2S7wRiM+wcPxL6T3PvNsM0VL993We5EFAAAYvCIijJOi00d4tn9sN46A\nx41D4Ea3U83SzqOe2/9FRUjfuY1j2wEACEfRkcY2v+Y4z/aDp6WDp/r+HniPwA1Jkv2S9PYnxkrl\nLhEytv/r+c4XAACED1OMdP/txifaV/tHvTHNFANH4IYcl6XKj6X2Hnttz75VmpDS9/cAAIDwkWQy\nQndMlGf7e59Jn30RnJrCCYF7iLv4pRG2L1/xbL/rZilnVHBqAgAAgZc2XLp3kjGdtItb0t+PSieb\ng1ZWWCBwD2EXv5TePNL7YJv8dGkqe20DADDkjE2QFmQb00q7dLqlt21SA6HbZwTuIaorbPc8Ver2\nVGnGOPbaBgBgqBqfLBWO92zrcBkH4xC6fUPgHoKuFbZv/fcvGGEbAIChzZJm3IC7GqHbdwTuIebr\nwnbRRCmSsA0AAGRML51xs2cbods3BO4hpOU6YZtTJAEAwNWmZlz7Tne9Izg1DUZErCHiglOqsBK2\nAQCAd6aOuXbo/vR8cGoabDisewg41SK9Y5O+7LHPNmEbAAD0R9fuZf/411dtLre081PJeUWakh6c\nugYLAneYO/6FcVz71SdISsaBNvMnELYBAED/TB1jrPXa+7ln+/9+boTub2Wy8cK1ELjD2OFGafdx\nY9P6q00eJd19CwskAQCAd/LSpfho6d3PjDvcXWpOSc52afZ48kVfCNxhyOU2PvL58Ezva9/KlKZl\n8A4UAAD45rZUyRQt7TgqXXF91X7EbmzQUJQtmWKCV18oYkJBmLl8xTiqvWfYjpBUeKt051jCNgAA\nGJhxSdJCi3G3+2oNLVK5VbJfCk5doYrAHUbsl4wXec+9MaMipP+cJOWMCk5dAAAg/IweIS3OlUbE\neba3fClV1ElH2cGkG4E7TBw9b7y4W770bB8eK30319iRBAAA4EZKMkkP5kpjRni2d7ikv39qTHF1\n9VxMNgQRuAe59g7pf44ZL+oOl+e19BHSQ5OlUcODUxsAAAh/w2KN6SWT+/gk/eBpaWud5Lgc+LpC\nCYF7EDvdIm2qlT629742eZRUbDF+CQAAAPwpKlKadauxXqznLiVnL0qba6W6Rsk9RO92s0vJINTp\nkqobpIOnem/5Fxkh3XOLlMt8bQAAEGA5o6TkYdKOT6RLV75qv+Iytir+vEkqHC8NG2K7mHCHe5A5\n2yptsUr/7CNsJ5mkBycTtgEAQPCkj5BKpki3JPW+dqJJ2vCh8en8ULrbzR3uQeLyFen/6qXD5/q+\nPmW0dNfNUnRUYOsCAADoaViM9J1JRm7Z97nnft1t/15/9vE541P5lGFBKzNgCNwhrqNT+uiscYJT\ne2fv6zfFSHMnSDcnBr42AACAa4mIMD51H2uWdh2TGi96Xj/VIm36SLKkSdMzpZvCeN0ZgTtEdXQa\n7wr/ecpzDtTVJo2U7s6S4ofYPCgAADB4JJqM/br/eUo60CB1XjWVxC0j73xy3gjnU8eEZ/AmcIcY\n5xWp7qxkbTT+vS+J8cZK4MyEwNYGAADgi8gI47TriSnSnhNSfY9D+jpcxinZ1rPG0fFT0sNrqgmB\nOwS43NJJh3TELh3/4tobxMdFGS/Wb4w2tt8BAAAYTBJN0v23G4sn//EvydHmeb3TbdzxPnzOOEzH\nkiaNT5FiB/katYAG7jNnzui3v/2tDh06pPj4eM2ZM0dPPPGEYmN7f3awc+dOrV27VvX19crMzFRp\naanmz58vSXK73SorK9O2bdvkcDhksVi0cuVKTZw4MZBPZ0AuXzHmLn3eZLzo2jqu3Tc60thX+44M\nycT0EQAAMIhFRBgnYI9LNIJ1TUPf02dPtxqP3celzETp1iTj031zfOBrHqiABu7S0lJlZ2erqqpK\nra2tKi0t1QsvvKDHH3/co5/NZtOKFSu0Zs0azZw5U/v27dPy5ctVUVGh7Oxsbdy4UW+++aZeeeUV\nZWZmat26dVq6dKl27NihuLi4QD6la3K5jY9HOjqNF9HFduOUpQtO6dwl45/XExMl5aYZ85k4wAYA\nAISTqEjjU3tLmnS4UTp0Wmpt792v023coPy8yfh6RKw0aoQx5STZJA2PNR6xUcbPDMVZAAEL3Far\nVUeOHNH69etlNptlNpu1dOlS/eY3v9Fjjz2myMiv/u+Ul5drxowZmjt3riRpzpw5Kigo0BtvvKFf\n//rX2rRpk77//e9r0qRJkqSf/OQn2rBhg/bu3dv9PcFyqd04Zv1Ui+8/IyFO+ka6ZEmVYpn0AwAA\nwlh0pDFne/JoY2rth2ekM63X7t/aLrVekI5d6Pt6kkmafYuUEUJr3QL2HuDw4cNKT09XcnJyd1tO\nTo6am5tVX1/fq29OTo5Hm8VikdVqVVtbm44dOyaLxdJ9LSYmRtnZ2bJarf59Ev3wQYNvYTsuSro9\nVVqUI/1XvpSXTtgGAABDR2SENCHF2NHk4TxpWoZk9mHiQtNlafeJG1/fQAQs0jkcDpnNZo+2hATj\nrUdTU5OysrKu27epqUnNzc1yu93d39vz+mAycpg0NkHKSjIWBoTiRyAAAACBlmSSCm6WvpUpnXca\n690amo0TtzsH4QmVAb2H6vbiDM/r9fXmZwXSN8dKXzilsxeNd2rRkcZCx675RSnDjEfacCmeO9gA\nAADXFBEhpd5kPL45Vup0SfZLRgi/4JSa24zpvJeuGOvmOlzGTiizbwl25Z4CFvmSk5PlcDg82rq+\nTklJ8WhPSkrqs29KSooSExMVGRnZ5/Xs7Gw/VO6dm2Kl7+ZKbrfxIgEAAMCNERUpjR5hPPoSqvkr\nYJMYcnNz1djYKLvd3t1WW1urlJQUZWZm9upbV1fn0Wa1WjVlyhTFxcVp4sSJHvO129vbZbPZlJeX\n598n4YVQHGwAAIBwFqr5K2CB22KxKC8vT6tXr1Zra6tOnjyptWvXasmSJYqIiFBRUZGqq6slSSUl\nJaqurlZVVZXa29u1Y8cO1dTUqKSkRJK0ZMkS/fWvf9XRo0fldDq1Zs0apaWlacaMGYF6OgAAAEC/\nBHQW8fPPP69Vq1Zp5syZio+PV3FxsZYtWyZJOnHihJxOY3PqCRMmaM2aNXruuee0fPlyZWVlqays\nTOPGjZMkPfTQQzp//rweeeQRtbS0aOrUqXrllVcUE8OpMAAAAAgtEe5QXX14AzQ0NGjOnDl69913\nNXbs2GCXAwAAgDB0vczJRnQAAACAHxG4AQAAAD8icAMAAAB+ROAGAAAA/IjADQAAAPgRgRsAAADw\nIwI3AAAA4EcEbgAAAMCPCNwAAACAHwX0aPdA6+zslCSdPXs2yJUAAAAgXHVlza7s2VNYB2673S5J\nWrJkSZArAQAAQLiz2+0aN25cr/YIt9vtDkI9AdHW1qa6ujqlpqYqKioq2OUAAAAgDHV2dsputys3\nN1fx8fG9rod14AYAAACCjUWTAAAAgB8RuAEAAAA/InD76MyZM1q2bJmmT5+ue+65R7/73e/U3t7e\nZ9+dO3fqgQceUH5+vu6//37t2rUrwNXCV96Mc1VVlRYuXKj8/HzNmzdPr776aoCrha+8Geculy5d\n0qxZs/TEE08EqEoMhDdjfP78ef3sZz9Tfn6+pk+frt///vfXfT0gNHgzzhs2bNCCBQuUl5enefPm\n6eWXXxazbAeHTz75RPfee68KCwu/tl8o5S8Ct49KS0uVlJSkqqoqbdy4UYcOHdILL7zQq5/NZtOK\nFSv005/+VPv379fPf/5zPf744zp69GgQqoa3+jvOtbW1euyxx7Rs2TIdOHBAf/jDH/Tiiy9q586d\nQaga3urvOF+trKxMFy9eDFCFGKj+jrHb7VZpaakSExO1Z88eVVRUyGaz6f333w980fBaf8f5/fff\n1+rVq/XMM8/o4MGDKisr05///GdVVFQEoWp4Y/v27frRj37U504gVwu1/EXg9oHVatWRI0e0YsUK\nmc1mZWRkaOnSpSovL5fL5fLoW15erhkzZmju3LmKi4vTnDlzVFBQoDfeeCNI1aO/vBlnh8OhpUuX\nqqioSNHR0Zo2bZruuOMO1dTUBKl69Jc349zFZrPp7bff1qJFiwJcLXzhzRjX1NTo+PHjevLJJ2U2\nm5WZmakNGzZo/vz5Qaoe/eXNONfW1mrixInKz89XZGSkbrvtNuXl5clmswWpevSX0+nUli1bVFBQ\n8LX9Qi1/Ebh9cPjwYaWnpys5Obm7LScnR83Nzaqvr+/VNycnx6PNYrHIarUGpFb4zptxvvvuu1Va\nWtr9tdvtVmNjo9LS0gJWL3zjzThLxtg+9dRT+uUvf6kRI0YEslT4yJsxrqmpUXZ2tl566SUVFBRo\n1qxZevHFF6/55guhw9s/s48dO6b9+/ero6NDNptNtbW1mj17dqDLhpcWL16sMWPGXLdfqOUvArcP\nHA6HzGazR1tCQoIkqampqV99e/ZD6PFmnHtat26dHA6HHnzwQb/VhxvD23HesmWLYmJiVFxcHJD6\nMHDejPHZs2dltVplMpn03nvv6Y9//KP+8pe/aOvWrQGrF77xZpzz8vL05JNP6oc//KFyc3O1cOFC\nPfzww/r2t78dsHrhX6GWvwjcPvJmYQWLMAYvX8bupZde0muvvaZ169YpMTHRD1XhRuvvOF+4cEFl\nZWV66qmn/FsQbrj+jrHb7dbw4cP16KOPymQyafr06XrggQf0zjvv+LlC3Aj9Hef9+/fr2Wef1auv\nvqqPPvpIr7/+ul5//XVt377dzxUikEIpfxG4fZCcnCyHw+HR1vV1SkqKR3tSUlKffXv2Q+jxZpwl\n4xd75cqVeuutt7Rx40ZZLJaA1ImB8Wacn3nmGS1evFjjx48PWH0YOG/GODU1tfuuaJeMjAydO3fO\nv0ViwLwZ502bNqmwsFAFBQWKi4vTtGnTdN999+mtt94KWL3wr1DLX9FB+a8Ocrm5uWpsbJTdbldq\naqokYwFGSkqKMjMze/Wtq6vzaLNarZoyZUrA6oVvvBlnyQhjH374oTZv3qyRI0cGulz4yJtx3rZt\nmxISErR582ZJUltbm1wul3bv3q3q6uqA147+8WaMJ0yYoIaGBrW2tnbP0W9oaOjXnFEElzfj7HK5\nes3L7+zsDFit8L9Qy1/c4faBxWJRXl6eVq9erdbWVp08eVJr167VkiVLFBERoaKiou6/fEtKSlRd\nXa2qqiq1t7drx44dqqmpUUlJSZCfBa7Hm3E+ePCgKioqtH79esL2IOPNOO/Zs0fbtm1TZWWlKisr\nVVJSosLCQlVWVgb5WeDreDPGhYWFGjlypJ5++mldvHhRhw4dUmVlpRYvXhzkZ4Hr8Xacd+3apQMH\nDqijo0NWq1Xbt2/XvHnzgvwsMBChnL+4w+2j559/XqtWrdLMmTMVHx+v4uJiLVu2TJJ04sQJOZ1O\nScbdkjVr1ui5557T8uXLlZWVpbKysuvuH4nQ0N9xrqiokNPp7PWH9Z133qnXXnst4HXDO/0d59Gj\nR3t83/Dhw2UymXq1I/T0d4zj4uK0fv16rVq1SnfddZfMZrN+8YtfqKioKJjlo5/6O87FxcVqaWnR\nypUru3eU+sEPfqDvfe97wSwf/bBgwQKdPn1aLpdLHR0dmjx5siTjkJtQzl8R7lCaUQ4AAACEGaaU\nAAAAAH5E4AYAAAD8iMANAAAA+BGBGwAAAPAjAjcAAADgRwRuAAAAwI8I3AAAAIAfEbgBAAAAPyJw\nAwAAAH5E4AYASJL+9re/afLkyd2P3NxcTZo0SR988EGwSwOAQY2j3QEAfXr22We1b98+lZeXKzY2\nNtjlAMCgFR3sAgAAoWfv3r3atGmTtm7dStgGgAFiSgkAwIPdbtevfvUrrVq1SllZWcEuBwAGPaaU\nAAC6uVwuPfLIIxozZoyefvrpYJcDAGGBKSUAgG4vv/yy7Ha71q5dG+xSACBsELgBAJKkmpoa/elP\nf9LmzZtlMpmCXQ4AhA0CNwBAkrR161Y5nU4tWrTIo/3HP/6xHn300SBVBQCDH3O4AQAAAD9ilxIA\nAADAjwjcAAAAgB8RuAEAAAA/InADAAAAfkTgBgAAAPyIwA0AAAD4EYEbAAAA8CMCNwAAAOBHBG4A\nAADAj/4f4JyXmYpWl98AAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import scipy.stats as scs\n", "import matplotlib\n", "import matplotlib.pyplot as plt\n", "\n", "np.random.seed(733)\n", "\n", "matplotlib.rcParams.update({'font.size': 14})\n", "plt.style.use(['seaborn-pastel', 'seaborn-white'])\n", "\n", "z = np.linspace(0, 1, 250)\n", "\n", "# We use Beta[3, 3] as our prior:\n", "prior_a, prior_b = 3, 3\n", "\n", "p_z = scs.beta(prior_a, prior_b).pdf(z)\n", "\n", "plt.figure(figsize=(12, 7))\n", "plt.plot(z, p_z, linewidth=4.)\n", "plt.xlabel('z')\n", "plt.ylabel('p(z)')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "z = 0.272052135336\n", "x = [0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0]\n" ] } ], "source": [ "true_z = scs.beta.rvs(prior_a, prior_b)\n", "\n", "print 'z =', true_z\n", "\n", "x = scs.bernoulli.rvs(true_z, size=30)\n", "\n", "print 'x =', x" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import scipy.special as scf\n", "\n", "# See https://en.wikipedia.org/wiki/Beta_distribution#Quantities_of_information_.28entropy.29\n", "def KL(a2, b2, a1, b1):\n", " \"\"\"Returns the Kullback-Leibler divergence between Beta[a2, b2] and Beta[a1, b1]\"\"\"\n", " return (np.log(scf.beta(a1, b1) / scf.beta(a2, b2)) + (a2 - a1) * scf.psi(a2) + (b2 - b1) * scf.psi(b2) +\n", " (a1 - a2 + b1 - b2) * scf.psi(a2 + b2))\n", "\n", "#print \"KL(Beta[3, 3] || Beta[1, 1]) = \", KL(3, 3, 1, 1)\n", "#print \"KL(Beta[1, 1] || Beta[3, 3]) = \", KL(1, 1, 3, 3)\n", "#print \"KL(Beta[3,.5] || Beta[.5,3]) = \", KL(3, 0.5, 0.5, 3)\n", "#print \"KL(Beta[.5,3] || Beta[3,.5]) = \", KL(0.5, 3, 3, 0.5)\n", "\n", "# See https://en.wikipedia.org/wiki/Beta_distribution#Geometric_mean\n", "def E_log_p_x_z(x, a2, b2):\n", " \"\"\"Returns the expected value of log p(x | z) over z ~ Beta[a2, b2]\"\"\"\n", " return ((x == 1).sum() * (scf.psi(a2) - scf.psi(a2 + b2)) +\n", " (x == 0).sum() * (scf.psi(b2) - scf.psi(a2 + b2)))" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimization terminated successfully.\n" ] } ], "source": [ "import scipy.optimize as sco\n", "\n", "def objective(a2b2):\n", " return (KL(a2b2[0], a2b2[1], prior_a, prior_b) - E_log_p_x_z(x, a2b2[0], a2b2[1]))\n", "\n", "q_trace = [p_z]\n", "\n", "res = sco.minimize(objective, [prior_a, prior_b],\n", " callback=lambda a2b2: q_trace.append(scs.beta(a2b2[0], a2b2[1]).pdf(z)))\n", "\n", "print res.message" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [], "source": [ "\"\"\"\n", "Derived from\n", "\n", "Matplotlib Animation Example\n", "\n", "author: Jake Vanderplas\n", "email: vanderplas@astro.washington.edu\n", "website: http://jakevdp.github.com\n", "license: BSD\n", "Please feel free to use and modify this, but keep the above information. Thanks!\n", "\"\"\"\n", "\n", "from matplotlib import animation\n", "\n", "def generate_animation(qs):\n", " p_z_x = scs.beta(prior_a + (x == 1).sum(), prior_b + (x == 0).sum()).pdf(z)\n", "\n", " # First set up the figure, the axis, and the plot element we want to animate\n", " fig = plt.figure(figsize=(12, 7))\n", " ax = plt.axes(xlim=(0, 1), ylim=(0, 1.1 * p_z_x.max()))\n", " ax.set_xlabel('z')\n", " ax.plot(z, p_z, linewidth=4., label='p(z)')\n", " ax.plot(z, p_z_x, linewidth=4., linestyle='--', dashes=(5, 1), label='p(z | x)')\n", " approx, = ax.plot([], [], linewidth=2., label='q(z)')\n", " ax.legend(loc=2 if (x == 1).sum() >= (x == 0).sum() else 1)\n", "\n", " # initialization function: plot the background of each frame\n", " def init():\n", " approx.set_data([], [])\n", " return approx,\n", "\n", " # animation function. This is called sequentially\n", " def animate(i):\n", " if i >= len(qs):\n", " i = len(qs) - 1\n", " approx.set_data(z, qs[i])\n", " return approx,\n", "\n", " interval = 1000 * 15 / len(qs)\n", " \n", " # call the animator. blit=True means only re-draw the parts that have changed.\n", " anim = animation.FuncAnimation(fig, animate, init_func=init,\n", " frames=(len(qs) + (5000 / interval)), interval=interval, blit=True)\n", "\n", " plt.close()\n", " \n", " return anim.to_html5_video()" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.core.display import HTML\n", "\n", "html = generate_animation(q_trace)\n", "\n", "#with open(\"variational_coin_closedform.html\", \"w\") as f:\n", "# f.write(html)\n", "\n", "HTML(html)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimization terminated successfully.\n" ] } ], "source": [ "import scipy.optimize as sco\n", "\n", "p = scs.beta(prior_a, prior_b)\n", "\n", "eps = np.finfo(np.float32).eps\n", "\n", "def estimate(a2b2):\n", " \"\"\"Returns an estimate of the objective based on Monte Carlo integration\"\"\"\n", " \n", " if a2b2[0] < 0.0 or a2b2[1] < 0.0:\n", " return np.inf\n", " \n", " q = scs.beta(a2b2[0], a2b2[1])\n", " \n", " N = 10000\n", " \n", " s_z = q.rvs(size=N)\n", "\n", " s_p_z = p.pdf(s_z) + eps\n", " s_q_z = q.pdf(s_z) + eps\n", "\n", " return (1. / N) * (np.log(s_q_z / s_p_z) -\n", " (x == 1).sum() * np.log(s_z) -\n", " (x == 0).sum() * np.log(1 - s_z)).sum()\n", "\n", "q_trace = [p_z]\n", "\n", "res = sco.minimize(estimate, [prior_a, prior_b], method='Powell',\n", " callback=lambda a2b2: q_trace.append(scs.beta(a2b2[0], a2b2[1]).pdf(z)))\n", "\n", "print res.message" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.core.display import HTML\n", "\n", "html = generate_animation(q_trace)\n", "\n", "#with open(\"variational_coin_mcint.html\", \"w\") as f:\n", "# f.write(html)\n", "\n", "HTML(html)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import theano as th\n", "import theano.tensor as T\n", "from theano.tensor.shared_randomstreams import RandomStreams\n", "\n", "from lasagne.updates import rmsprop\n", "\n", "def log_gaussian_logsigma(x, mu, logsigma):\n", " return -0.5 * np.log(2 * np.pi) - logsigma - (x - mu) ** 2 / (2. * T.exp(2. * logsigma))\n", "\n", "def log_beta(x, a, b):\n", " return (a - 1) * T.log(x) + (b - 1) * T.log(1 - x) - (T.gammaln(a) + T.gammaln(b) - T.gammaln(a + b))\n", "\n", "q_mu = th.shared(value=(float(prior_a) / (prior_a + prior_b)), name='q_mu')\n", "\n", "q_logsigma = th.shared(value=0.5 * np.log(float(prior_a * prior_b) /\n", " ((prior_a + prior_b)**2 * (prior_a + prior_b + 1))),\n", " name='q_logsigma')\n", "\n", "rs = RandomStreams(seed=733)\n", "\n", "objective = 0.0\n", "n_samples = 10\n", "for sample in range(n_samples):\n", " eps = rs.normal([1], avg=0., std=1.)\n", "\n", " sample_z = (q_mu + T.exp(q_logsigma) * eps).clip(0 + 1e-6, 1 - 1e-6)\n", "\n", " sample_log_q_z = log_gaussian_logsigma(sample_z, q_mu, q_logsigma)\n", "\n", " sample_log_p_z = log_beta(sample_z, prior_a, prior_b)\n", "\n", " sample_log_p_x_z = (x == 1).sum() * T.log(sample_z) + (x == 0).sum() * T.log(1 - sample_z)\n", "\n", " objective += (sample_log_q_z - sample_log_p_z - sample_log_p_x_z).sum() / n_samples\n", "\n", "updates = rmsprop(objective, [q_mu, q_logsigma], learning_rate=0.02)\n", "\n", "optimize = th.function([], objective, updates=updates)\n", "\n", "q_trace = [scs.norm(q_mu.get_value(), np.exp(q_logsigma.get_value())).pdf(z)]\n", "\n", "for epoch in range(250):\n", " optimize()\n", " \n", " #print 'Epoch {0}: q ~ N[{1}, {2}]'.format(epoch, q_mu.get_value(), np.exp(q_logsigma.get_value()))\n", " \n", " q_trace.append(scs.norm(q_mu.get_value(), np.exp(q_logsigma.get_value())).pdf(z))" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.core.display import HTML\n", "\n", "#html = generate_animation(q_trace)\n", "\n", "#with open(\"variational_coin_sgd.html\", \"w\") as f:\n", "# f.write(html)\n", "\n", "HTML(html)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.12" } }, "nbformat": 4, "nbformat_minor": 1 }