{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\LaTeX \\text{ commands here}\n",
"\\newcommand{\\R}{\\mathbb{R}}\n",
"\\newcommand{\\im}{\\text{im}\\,}\n",
"\\newcommand{\\norm}[1]{||#1||}\n",
"\\newcommand{\\inner}[1]{\\langle #1 \\rangle}\n",
"\\newcommand{\\span}{\\mathrm{span}}\n",
"\\newcommand{\\proj}{\\mathrm{proj}}\n",
"\\newcommand{\\OPT}{\\mathrm{OPT}}\n",
"\\newcommand{\\vx}{\\vec{x}}\n",
"\\newcommand{\\I}{\\mathbb{I}}\n",
"$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"\n",
"
\n",
"\n",
"**Georgia Tech, CS 4540**\n",
"\n",
"# Lecture 24: Sampling and Markov Chains\n",
"\n",
"Jacob Abernethy\n",
"\n",
"*Date: Thursday, November 21, 2019*"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## A Core CS Problem: Sampling from a Distrbution\n",
"\n",
"There are many problems on which we need to be able to generate random samples from various probability distributions. \n",
"\n",
"- **Cryptography**: Generate random prime numbers\n",
"- **Optimization**: Optimizing non-convex functions using *simulated annealing*\n",
"- **Physics**: Sample states of a physical system, i.e. the spin direction of particles in a solid\n",
"- **Machine Learning**: Sample parameters from the posterior distribution of a Bayesian model"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Sampling has had a rebirth in ML\n",
"\n",
"We have recently seen the birth of Generative Adversarial Networks, a deep learning technique for producing random samples that appear to resemble images that match the distribution of some training set."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/jpeg": "/9j/4AAQSkZJRgABAQAAAQABAAD/2wCEABALDBoYFhsaGRodHRsfIiUlIiIfHyUlJSUlLycxMC0nLS01PVBCNThLOS0tRWFFS1NWW1xbMkFlbWRYbFBZW1cBERISGRYZLRsbLVc2LTZXV1dXV1dXV1dXV1dXV1dXV1dXV1dXV1dXV1dXV1dXV11XV1dXXV1XV1dXV1dXV1dXV//AABEIAWgB4AMBIgACEQEDEQH/xAAbAAABBQEBAAAAAAAAAAAAAAAAAQIDBAUGB//EAEIQAAIBAgMFBAcFBwMEAwEAAAABAgMRBCExBRJBUWEicZHRBhMyQoGhsRRSksHwFiNDU3LS4RVi8VSCk7IHM6I0/8QAGQEBAQEBAQEAAAAAAAAAAAAAAAECAwQF/8QAIhEBAQACAgMBAQEAAwAAAAAAAAECERIxAyFRQRNhBCJC/9oADAMBAAIRAxEAPwDz8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAC9HZU370PF+Q+exqqtdwtzu7fQzyjXDL4zgNin6OVpaTpeMvIf+y9f79L8Uv7Rzx+rwy+MQDZfo3WXv0vGXkQV9iVKavKdPxl5DnPpwy+M0CWph3HWxG4l2zcbCASug+aD1D6DcONRATfZpdBrovoNxeNRgSepfQl+wztfL5jcONVgJ44WT4r5hLCSWrQ5Q41ABNDDOUlFNXbS48XY6Ct6D4qEbupQ5JKU7t8l2RuJquZA08dsKrQdpypt/7W+duKIqWyqkle8V3t5/InKLxqiBehsmpJ2Tj4vy/wADa2zpQ1lB87N+Q5Q41TAe6bXFDDTIAVRNfCejlatFSjKmk1fNy8iWydrJb0xwOg/ZDEffo/il/aRv0Wrp5zpL4y/tJyi8MvjDA34+iNd6VKH4pf2jZ+iteOtSj+KX9o5Q4ZfGEBp1dh1I6zp/By8ipPByjq4+L8i8ocargPlTa1aGFZACpCWAAHqnmldK/F3sOqUXF2un1WaC6RALuhug0QCWOHb4omhs+b4x8X5GblIswt/FQCy8FLnH5+Rb/wBBq/ep+MvIc41PFneoywNyh6LV6iup0vjKX9pN+xuJ+/R/FP8AtNS7Yssuq50Do16GYn79H8U/7RV6FYr+ZQ/FP+0I5sDpV6D4r+ZQ/FP+0ZW9DMVBXc6L7pS/tA50B9Wm4ScXqnZjAAALuH2ZOok4yhnzb8ib0slqkBsQ9HKz0lT8ZeRIvRau/fpfil/aTlGv55fGGB0C9EMQ/wCJR/FP+0R+iWIXv0fxS/tHPH6fzy+MADbl6L1179L8Uv7SGWwKq1nT8ZeQ5ROGXxNhKl93o7GxDQwcO8+Wea8jcoyyPPk9WJXGUPY8Hp8B8cZK2cWKxplolSvJ6KxSqQbe9J3ZamQziBz+Nleo1yItGh1bOcn1YzieqPNfqxFaeI6KGx0FTzMtJYuxHVjxWjJsLFSk0+RFVybRGoZHUvwVoXksnoyrh57ru7fEtYrFRmsslxXXoZreKtBW+IzE1LiXzIqnIs7Zy6LhHetT/rj/AOyPWNoVLNv7t0u9/wCPqeUYSP76l/XH/wBkeg7cxu7PstPV34t65fW5c+nLDth7RSVRznnraPdll5/8mZW2jk1GMbX4L9XG4uq5K7ecs/hol3FKTWi4GccW8stJJ4ybVnLdi/djZfQq1arYyTBI6yac7ls3gIotuw7d1Qscs+PDzKyWpbRaL5vmdvsGd6Ue5HESjkdV6NV7wSMZ9Onj7dKRVYXRKtAkji7MTF05Rzi2u5mXWxtZZN3+B0eJhcyMVRRpGJXxk+KRRqVpPoXsVDMz6iOmOnLPaJjR7GHRypRWNAIdwFhUa7hogXZ7fIIjR8WrEqxdpRTQ/wBZwWhDRnawreZws9u8vpPKWTOhisl3I5WtPJnYulaKvyX0M5enp8OW9n4Gruz6M2EzAeRsYOrvRR08WX45f8rD/wBRZRIiOJIjs8SWI6pT3kNiTwQHlnpPg/VYueWUu0viY52v/wAgU4J0/vt//k4uwAW8Bi3Tln7L+RVBonay6drhK6a1LyqI4XC46pTyTuuTNCG3JLWPzONwv49OPkddGt1EniElqco9uSekfmSwWKxHsrcjzeRm4/W+W+mhjtpxjlfPktSKhhqtftTvTpvh7zJsBseFLtTe/Pm9Ea8FzJvXRr64VO3a4o3cNnFNaBU9H4WSVSd+5FjDYFUoqKblbi7KwysZlLYa0PlUiiJ11y+Zldwu6RYpbsG+g77S+EF8W/IhxVT1kGpWgnlk7lOUc5vZdbiOJdeBis1J26pIPsafHI77jjqqcJkkmS/YVfV+Ba/0xP3nbuJbFm1OnPPLkIld/rQufYlF2u8+NiKrDc6k2vRteitxSTzvYrwvfoPliG1ZpeIkcRp2fma1U5S1JK3DIZVWaF3r8B8Y31JqtXKUlFfvKf8AXH6o3ttYnfr1Kkk7exHPkZlHCq6lvPJp6csx+KxPrHeyTu2viTVyZ5SKc5JJRzcuP+1cl1IoQzJVRSd7scoWvnr8jeq58paoVFmDLc8Om73Y37KubNe2dxE43jcSnC7SLMaFk1cI0Ene5F3FapG7b6mhsHE7lWzeT+pDKim9RKdHdaaejuSzcWXV27/D1lKKJmcjh9tTppLcTt1sWV6US/kx/G/I58MnX+mLeqxuZOLpleXpPJ/wY/jfkVau3HL+Gvxf4LwyP6Yq2LpmZWiaVbGb/upfEpyimaxxrGWeNUZIYXnQXMa8MubOmq5bimKkWvsq5sX7KuYFViIuPDJq1xPsa5sCoBb+yL7z8A+yLmwK8J2JHMk+yLmxywq5szY1KjnbdPRNqRS3LcYr6Hn86C3XnwPVK+ylWjBubjaK4X4HPPHeOnXxZzHPd6cxIubNnm0ab9HI/wA2X4UPoej6hK/rZP8A7Uc8MMsa9Pl83jywslJAkii1HAJe8/Ad9kS4/I9O3z9IYlfam1aeEourPPhGN85S5IuypWOZ9Idi/aqiqSqyioxtGKirLn4k5RdVxO1NoVMVWlVqPN6JaRXBIpM3cRsKMP4kn8EUZ4CK99+CLuLxqgKXVgI/efggeCj95+CG00mwmzIzV238DRpbAp2vvyfgVcNivVRsop/GxYhtmS/hr8T8jlly/HoxuEntfwuy6VOSajn1zNCUjBntuS/hr8T8hv7QS/lR/E/IxccnT+mDfTJYs51ekMl/Cj+N+QzE7cnVhuqChfVqTbsThU/pi6iraEdcvmZOIr72WkeS4lnF1t7u+r8ii1fueWWr6IxGdI5Sf64D4U7ZvJfP/Asmo9/BLRfrmQ1Xxk/gjaUs6uqivyX+SnUqZ5dp83ou4KlZvsxXh+sx0YKGc3d/I1IiOMcruzZHVq9Ba+IvlEqykakZtXKNWL17i/Ta3ctDDpyd8rXNjBRcumRMo1jU1SlZZZozcZDI1oylB2a7P5FXa2H7O9Hjo/yZnG+zLpz0mMuPqfMjZ6HBcozViaLVzOjKxJGrYzY1K2aclYrFRYrRFplwmkzuyAAHRzAAAAKIKZqgS4NlerX4IipZTS1GetXUr34hF3kNmk7rdBYTvyIJZMmpRV89F83yG10tQjF+98mLKlH7y8GJTr5aqEOC4snk7xuvglfPvJumlKbS4oZ62POxPUp3V559EVZq/JIsyqaSqVxxSTa0ZPTq311NbTSYVDbi3LpJTrBYaKjOl2dYBA3iNEqaPuZ6/R9iPcvoeOVamT7j2SguxH+lfQgfYWwAUAlhQII5xMnHxsmbEjntt41QjKT7kubMVvFy+1a9pNIyVzY+tU35NvMbe3ebkW0rYMIR4vwEqyCGSGxQ0npxyAimQWJKssyO5A9LIEhY594nEDq503LtS0vZLmyOatnlvP5LyLtbKzaz0jHkinVainKXfc8sehUqzUFvMoSlKo8+zH5vvJ5xdR7zWXur6skdOytxOkYs2qzq7uUUVm3LXPoa1DZrnrp08y8tkqKyt4DlF4ubdCXJkXq2dDW2dbR5cuZTnQkldKxqZM3FkqNmbOzamVjOnCW9lc0cBGS1fiTKkjWxdBOO8u8zZVV7EtHkjVozTjbQxNpQ3W1w/VmYnbd6Y+Po7snbmUzSry342eq1ZmyVmenHp5siAAGmSx1XeajMuOq7zUZYlIAAVAAAQDEAjqVLaEqm15ZWRXUSXVcx9KLTM7aNhD2b6O6BL1dRPh+TLfqlOF4acVyZBu7ySeqvb4EDo0leUJa3yIqUbz3X+syao7uEuln3kWd78bFCVau/Uv7qyiuhNTne2ebyS5LmV1Bpsluoq2fw4gWHeVoxvbnclnQtkrOXQqetm0vdiuXEtRxO4t1R7X6zfMmhQq0uV/ArtW6GpPeu2+1zUb2RWq0k7y06IqIqVW+upOmUXkyxTnc3tLE9wlOwiEnG6CRFOvyIXUYSg0NsZaDZ7hQ9iP8ASvoeHtHt9D2I/wBK+hFiQUQUAAURkVVxlWytfU8+2/j/AFlVxi8lkjp/SbH+qoyfvSyR5/dt3evHp0Mz3WuofH6CxzYkskktWPhlpr+RsLN2XX6FeWZKs83oiOYCRRI3l5DVGyG1JdlmdrpFPNkfGw+Qk0VD4PMfJEUSZvK4HW1M22+WdvoihUvVln7K18u4vV08orO/6uOo4ZacFr1Z5NvTpUp0L52/x0LGGwd3mXoYe7LlOioi1ZFenQSQ9wLDiMaI1pWnTuilicOn3mlJEFSBNmnNYqluvQZTqOOauunA18ZQvexktOOp0l252L1Gto3+ug3aFHeg7ZtZ96KlOtaVno/0i5Qq37L1WnUt9e0c5V5lOoszX2nh9yV17MvkzJqqz/M9GN24ZRGwFYhtzLHVd5qMy46rvNVlKaAAEACiEUkpWVyspXeYuInwIkyUi9SpxtdW8QnOE8l2X18yClUvpZFqGG39DFak2ghvwlda/UsQpKorpWfFGnhdmOSzVkamH2RBcLmLm6Tx1gU8A5cMixS2XfJ5aNM6Wng0uBPDCRRjnXT+ccrHY4lbYz4HWuiloiCpSQ504RxssDuPi/gRuGefidXUw6Zm4vBZXSNTyM3xsSriLKyRSk73L1alYg30lmkdZXCxQlxFhKzJJxyutCJGkXou49Feg8u4sIqaI4ierXIcBmtI6lJbr7mexUPYj/SvoeQT9l9zPXaD7Me5fQlqpxREKihSOq8iQxvSTaCoYac79prdj3szSOI9J9o+uxMkn2IPdXWRmQX/ADzfFkSV33fVkr5Fi0jfi8kSLJfQiWby7ia63orgs/Iqlmt1JePUjjG7uSyV/iTwoWsZtdMYhVGyz1ZVrRtkblKhdrIzMZStNmJWssVG2Q5L6Em6CWh0c9IVHOxMlp1FnH9dUPUez8ybNOyjS7V9W/1cnp0+mQrWduL+hZpxPHHohYU7Eu6CHIqmboySJmhkgqvKJBNFmRDNEVUrQuY+Joam3NFPERuJdM2Ofqq2pNRqXWvaQuNpq5Xpzs09La+Z2nuOV9NKpGNak0+OT6PmcviYOLaeqdjo97d7Szi9UZG1aS3t5aPj5mvHdXTHknploAA9DgWOq7zVZlR1XeajBSCiCigGTlZDrkFdkFdyBZsJodSjxIq1hqWaOk2fg+LWRkbNob8lyOtoRSSOHkyejx4pqNBIsJEcWSRObto6Mbk0IDYksCpUdSBVqIuT1K1VEIqVEVZot1CvMjVjKxeFRgYmhaTXA6usrqzMLFqzdzrhXDyYsSXZYbt1dE1SGT+Q2Mew2d3mNw7zaLSZVpe0i0jRShYAJVhJ6PuZ63h32Y9y+h5HPR9zPWaMuzHuX0OeTUW0ORX32LGvZ2Y5Gk83ZHn3pzjd6vCin2YR3pd7/wAfU7utVW7e+SzPJNpYp1q1So/fk7d36sO6RFTeV33i72RHKXDlqEnwNixR/X67hm/x5sRysrfD4vUbHOSSCxbwsHJpsvPJjcJSsFN78ujdl3HKvRjNNLDRyTfEpbRo9p5ao1KMOXAi2hSur8jLenPVKdvkMqQs+hcqJPL9WIHrbkbjllBu6cn+rhGFrrwHxV7xWqzQ5LeVzWmduww13m9WXIIxqOOlHWJfobRg7XR5dV220EhUhITT0Houk2a0RTJZSsZeO2juO0Vnzeg0bWpkE5Iy3VrVdL26ZDlsurLNyt3tjivJZnNcypVeY6eypLNSz7ytvSi92efKX5MnFdquOhePMxHV3ZHQYiN0zmsZlJnTxuefr22sJUUo7r5ZGZjezePIbga7i10Ze2tR3rTXFJmpOOTF94sGceQxktVWZHJ8Tu4COq7zUZlx1XearKhAAREAyrUd5FmTKsn2hUIlrlxH6OxYoVo7zUrZqwyor1UuC0MrG7siFkb1JmNs32TYpWfE82Xb3YdLECxTRDBFiGRGksIksUMpse2isU2oitWiWnJEFeXILFKaK80TzeZVrYiC1kl8TOmtq9YyNoQyvzNSeIhLSSv3lDE55G8fVZy1Y5+crIilPs2LGPpbrvzKTkeme3is1Tqb7SL6jczU8y9B3NIksA3dFuQJPR9zPWcMuzHuX0PI6k8n3M9gwq7Ee5fQ55xuVLujJwTJrENbJN+BNDB9J8b6nB1Gn2p9iP5/mecX4vgjofTLH71eNGLuqaz/AKnqc3Ulw8TUhadF8X3jou2ZFBiuXE2h0pWL+z6PveBnQzaNnCYiMdU7cDGTph2u1OxSk+nzeX5jdkUt6fSK+Y3FYqFWO7B5u2uRpbHw+5Bt6s53p33ur9GlYWpRumuZNR0F4kb25XHYdwqOPPOPkUp5xvxTs/yOr2jg/Wxy9pZo5qrTcJN26NfkWVnKbitTq2+Y+U92V46OxHiIJdqPsv5Mip1OelzrK8+XbvvWU3rb6DVg433o5GFXwVek5L1jnGL3ZNWyeul9C3srGySW/o210y5deh58sLHXHOZOgoK2RZWhThNaosufZMyt1BiJlOOFUnd5vqJiq2Y6eI3KafQntKnhBRI6uJinbeXiZmKdV0/WtuNN5LLtPm7cNDFcak5uLm4ve3Xo1fhmtdTrMNsc5HTPFIrYqO8rmBCtOlVcJvR5NaNG7h5NwvzM3Hi1Mtqk45HN7Tp2mdZUjkYO1qfEYXVMpuMzCxzNqecN37quu55P5/mZ+zqObNCatUiucZL6ebN5X2xJ6YeMp2b8SmaOJztfV3T77eZnHadOF7LHVd5qMy46rvNRmkIAAQNmU3qW3xKsr3YqQSauybAq9RFWRb2cv3iM3pqdtuCnKThDJcWXfsMoxyrZ9UVZt0m2l7RnVsdKUrXk3okjlN3p3tk7aksXXpZb0u/KRNh/SGek43+RjOVVS3ZqMGrX9ZJq1yTDzbb3otLLP2o56X4otx+pMvjrcLtDfVy79oOdoJ05JNW+jXQ14Ps3ONeidH1sVu5vQxsX6QOLaivlcmxdXey1MzFUt1XfZXzfcaxZy/xXqY6vVfFLwJobL96rNLone/iV69CtTpqp6pQg3ZOWcu+xSxE6k5JbzbySW6o/Q7arz3X77a8tm0X7MnfvRDuyhLck95cH5mXQxMoT3XdNcMzZpxc1dmbLO2sdXpl7TXZ+JkG9tWHYZhuJ0w6c857NRYp1bIv4XYkpQ3nK3QzasHCcov3W0XlL0zcbjN1M6jtoNVW5EpjLlZSVGz2jDPsR/pX0PFN7I9oo+xH+lfQxk1FpyyMzHY3chUqPSnF+NixUm1ozk/TLHqFNUIvOecudr/mzM9tOPxFZzqSm8222Vb3FnIajrGac5WEuNuLwCL2y6O/O70R0TlBRSaVupj4CDjS78xlqlSe6vnojnZuu0vGNCph6U/ZaT6EsMNWhnGd1yMqrQdOTW88nbp3mrg4TdOUvaUHZyin49ULjSZe2vg8RKyUtS/TzMahX0Urd/BmzhFdI5u2z5RMPa0VKVt1X4s6KtGyMnE6t2MtOXxuGdPPg/wBZmep2v1N7ExnWjNxtuwTcm9F0XNmCqE3v5x7Nm7tK9+EVxO2PTz59uzxmBrV5b0txPK7imr9+YuH2ZUSjF1Mou+i8WbuSInmcbnXaYSdI4qzJnPIjcRKmhhvTPrO8glh5VFGKkkk7tNXvyEazLNNcjXTFNxVOtKnuXja6s7ZrqjLnsyre+/pnpodBGYkkXnThKwsNsq0t6Tcn1NNUkkTqA2SM3K1ZNKdaJz+0lvTUep0GJeRiVY70/j+vkTHsy6GFoKMb8X+SYzaD3cRS7vIvNdpRXup+NjP247Vab6fkjePuud6UMVG2/wD7bP8AL8kZLNXHT7U+sX/7f5MlnpxcMix1XeabMyOq7zTZayBGAMhTOZBVWZMyCpqWsn0MPv35IuYDD2qxi+ZDgZWTXU1KSiqkXF3s8znlXbGfrfqYVThZow5YGVCrvwSy5nR4eV0iWVFS4HGZWPRcduaxLeInGUoqMnZS4p9e80tnUXCLhGMXve1Jp+HcaMcHG+aJvVpaGudY4RR9Q4xtJp2eVlouRapL93YbVZJD2Dnt2/Gc4Z/ESphFOSnnloP95lqkiw0z69KcoODlJxeqeZmTwsou8ItPOztodOoDZU0b5VzuG3M4TZT33Opm+b1NKVNRRdqZFKvIzbutzHTI2krwaMvA0N6cVxvryNjEq43BYBwtJvNpOx13qOXHeSfDKXrN2TeTOe2nFKvUt9466MEnKo8rRz+P/DOLxE9+cpfebfiTx91PN1IiYgosUd3mJY9oovsR/pX0PH1DsvuPV/X2hG33V9DGTUOxVZQi5Ss1bQ8u2xjnXrzqPi8u7gdR6T7VtSdOLzllJ/VI4hsSLSNiXAEjbJ0EPpQ3ppdRryRNgn+8j3olWOkoUOyiGeEkpqUXZmnh45It08Pc4709PFkrAzqO8klpdrK/Wxq4XDyjT9Wpbq/2pX+Ny5ToErVlkXknGMSeE9Wrb182725m1s2PZVzOxKztzZrYRWSM1qJcRHIz8dh96naLs3q+hpyzI5U7mVYFClUpxcVuyi1ZprVFH/Qk3JtSd9EnZLM6adBchVTNzJLjKmcRGT7oySOFjqg0GtXQ9RWs2kursWKkEkTRtg4hWZYwk94MVKnvqm5JTeajfNhSpOlVjylkVF1RFsTboySCopENRk02VqpBRxjy7yhRV6t+Wf68PmW8VLPu+pWoLKUuj+eRY55Fpu8u+/5md6RS/extkl/g1MJG84/0/r6mN6Q//Zbqzph25ZdMvFyu4vovovIqy1Jq3sr4rw/5IZO7PTHGiOqNJmbHVGk2KkAjYEc5WBQ2QSY+4lWOhAUZWZ0OHwq9W5+9qc1c6/Bz/dQSt2oxbfwzMZ9Ovi+LODrZI06VQwcLKza5Nmth5nCvVj7i+hJjacsh8nkRdM+rLt2LsY9gy514xqXlzsasMTBQzGkrIqytIu4eRn4icZSaTLOz5XTT4FGlEjquwqeRDVmRVavMz68yzXkZ2ImXFcvURzz7uJcwkXJfr4EOFwzqPdXI0sKrQk37r3X3o3XGXXtmekGLVOiqUfanr3afr4nKSZb2ji3VrTnwbsu7gU2ztjjqPNnlyoHRGMVM0wtw0O6xON3YLnbwy1OGoyyNzaeL7Fr8MjLcY208Q6k+i0KLJJu7IzUSmj1khgFQrZJhZWnF9V9SJiwea7wbd3hHdI0qJjbMqXijZos872T3FlMjqMdcjrzUYNkFF9qokbWGjkY2BXbz1N+hSy1CEYWEqUxm+1qZU5iJdB6zHqIUiEmKJIiolFDql7XGjK1XKxJWtbVKtOMZ+scU5WsrirttN8CStPe7xtMlqr1N3Q2ZHTkOkyMoZleoyxIp4mWQVm1e1O3e2M0pSfPyDfzlzzXxtoFX2GuhXK9pdnvtL+n8kYfpE/3/AMZfU2Nn5JP4fXyRh+kLvXduj+SOnj7Yz6ZU811Tf6+RCyRvMjaPTHnpY6o0WZ0dUW61Xgih1SpYrylca2NuQSJklR5EFyS5Ay502y4uVGD3uzCLV+blJnMnSbEknhZe1dSa9qyzXiSzcaiehe7u03fVGnh5GJgp/vJxu273zNWg7M4Zz29Xiu40YyJosrQkTp5GHWs3GYGM5WkrxfyKtTep7sFdrRGtK183YcnB2zTKyxcPgFGfrJLen1Zo4Kk4XbepPOCu7NMY2TapXMhqTGykRykTSq9dlGcbst1syGSOmLGarW2nLDTg0k073Wl13kOI2rVqUq24koO0pvjmlF+L+pQ2pU3qzX3VYub7WzHF6OreNvvXzlL4Ky7jtJI8tyvTFYwemNNuYEACokhIvY+tdpX4GdFks5XZlZTGIxyBoojYCtBwKC2Vxo+WiQwI6bY9e8YnR4epkcdsKrrHk/qdRhpHDL1Xr8d3GnFkWKhvRaQRnYVzuYbZ7nKD0zRp0MW2klqRQjvMt+pSs7ZgMVTEb93uKnyzcvEtxi5K7Q+I5yDNqOGTsSOQySGbxGkpBiMXCGTeZLvEVWlGXDMy1jJv2pVMbfhK3RESxis/aXS1i2qe7qhst3LJdczck07emdVxNnlFvqsx8Ma46p27i3Jp6RVyOOFzuzNkLrS1RmpJNEkmNjkNnIw5EmzMxlTJ95brVDOqvelZBKr04Wt0z+LHONossugkkuL+hVryU+wsr5vohj7YFB7sId1/hd+Rj7Xpqc3KLu45SXHvNvd7LlbJL5LQ5utVfr3LhJ/LQ7YduebOmg3bx7ifEwcJWEwy3m480/p/g77cdKyJJMY1mOZWSCABQo5MaKiBbG1saeW7z1/IxkjU2RLtfAlWHuoqeKTTVnkzepM5PH39Y8819ToNnYr1lOL46PvOWc/XbxX8a1OQ3GbQjSjduw2kwrwi1nFPvRyehQoYx15tQ7TXPJF2OBrtXW5l1/wLSoU5XaSi+mRYhS3IOO+916ZvzN7jVmX5WfVwtaEd9rJcVJFWntKWfG3NZmtVw6cUpTbjyvkZ+Tk4wWQ3DVk91PQr760fxRLNWFp5IjrTMIhkVsVVUU2+CJpysjI2hUcoTtokbxjnndRkSlduT4u5fxErYaMU3bJ5638s2ZiZPVrOSiuSO9jybQsaOY0rIEFAoESMjJbEoEOULoko4aUtE7D50HF8uV+JNtaVd0SxJJETKhJPMaKFiizs6ru1Y9cjsMJPQ4ZOx1GysXvRT8Tn5J+uviy16dE1dKxTxONVN2eTLGHqXViPE0VK6ZxehXhtN8LFyntVNJOOfO5SpUt15pNdxajg6UlfdV+n67jTrJNLsNp21sSw2nTersU/sNGK9m/6/wCAWAhN9mCiudhYlxxXnj6d7b6u+pInfMiw+BhDRIsbtjFc/X4jpVVJEiZkbIqP1UG+SNWLMNJVG42WHjyQ6LHMoh3EtCOSJpRIJsyprIKtQWrUsZO0dpRpLnPgvMSbOk2KrWsrq7dld2uOwlJbzzyis3zfM57Z85VK0q03vOKy5JvJWXj4HQx7FFQXtTebfJGs5r05cuSCrWc5S3corV/RLqJQwuV3kr6cX5lijRT7r5dXzK+0NofZ7qKUqrWSayiubX5GZN+otqDa2NUI+pjbel7T4RXmYboprd4rPzRbVPfoSnK8pOTu+Ol7lKpU/eJq2dl8bWO+M055I8ZHsri1eN1xto/BkOzf/uh3r6k1epeM43vo1fXu8PoQ7Pkoz3n7qb8Mzp+Od7QVPafexrCTzBm3MjEFEKFHJCRiWYRW6QRI1tiU7t9F9eZkyWZuejML+t/7V4v/AAZvSxmYmF6tRcdV9S/sm6pNr7xQxck69R5rtNGvsmFqaM5X0345/wBmjh6t8y2jOlTcHdacUXsPVUkjhXqhtTDyTvDUP32XZTt1/wAGlStYmuixrlYxZUqs9bRXTMfGhu8DUnJWKlRkpvatIgqSHV6qRm1sQ5O0c2WRm5aJiark92OrHTwq9U480T4TC2zepYnDI1v4xrfuuOjhpb268htSm46/A6enhIznaSytrxS/PuMnbUouSSte704I6zPd08+WGmWxBWhDo5gQUAA7Wl6OQppPOUmlmzij2fD0FuRdvdX0JVnph4HYUVm73GbX2NSjhpylk43aZ1MIWOZ9NMbu0lSWrd3+RmzTUttee1lb4FcnqyuyBmolBJGGS7rjbZE2ke9WCKxZwGK9VO/uvXzK9hpSXTt8HibpO5oU3vM4jZu0HSe7L2fodZgMSpJO90cMsdPTjlto+ouKsIS0XctU0YdJdK9LD24FuMCaMUDCW7R2GzY6TKtetYgp7NVopF7dtoUtnLKxopHJsxSkI6rXAlsI0aNoZYro/Aq1sTyT8C3JEFWFzKsfFYmXBZnP4l5tyd3c6LH2jFs52rC7u3kdMGM17Y9JylBcG7/kbG9v1G+Gi7l+vkZmy52cpaWjkX6fZjfpx4mPJfbOM9JsfjVQVku3wXXy5mHCLcnUm7tu7fO9v8i4iq5yc5Pu7l/wNjUvFrp9cjWM1FOnPcpTiuDT8U1f5GFVl2viaMq11nrp5fMq+qUW3POX3eXe+Hdqd8PTln7Ryg5Se6m7p8OgycdxOL9p69ETOq3ley5LQZioWSfM1KxYqMdYbYdc25kYiAWKKJaKHO6GJkl0ZEUjc9GH+8nfRpJf1Z2MK5s4a3+nzlHKpGqnfi1kvMXpdqWKp/vZRtm5P6nQYKnaKRj4GMq0/WTzfM38PE5Z38dvHP1I4lSadN3Xs8ehoboyUDk7Fo41NLMn+0mLiMNKLvDwIVVq6F0cm866KmJxaXEzd+o9WySlh7u7uxo5fEc5TqvLJF3C4RR7ySnSsTpC0kKojJofcbYkWq0odE3ZrPnbJ34Z2MHGYNQhvSlvTm1Z3u+rOilkYOPi1Umnxk/hnodcHDyMmrTa7iI05Q8GRPDRfCx2286iKWng+TE+ysbFax7fh49iP9K+h41Khk+49noexH+lfQKe8keaelWNVbFSSdoR1t3X+f5o7T0g2j6mk4rOTTbt93/Lsvicj6ObMWIxM6tTOlTblNv3p8u6+fdYlI5zEUZRsmrOy4Wy4EDhZX8DpfSmhuV4JrNxc31bbaXysc9XlmFRNEzd4LoRapEkfyZUMcFut37iFocxrKENDZu0nRdndw+hQEHZLp6Ls7HQnFOMk0akKh5hgMQ6dRWk0m87Ox2dBVGlao/ikefPHT04ZcnQxrCyrLmYlP1qdnK5ajS+82zG201XEt5RIdxv2iWMB7iQQ4JfQ0UjlqXpJTjpTk/iiT9s6a/gT/FEk8eXwvkx+umsJJHM/txT/wCnqfiiJ+21J/wKn4omv55fE/pj9dFIimjn36bUv5FT8USOfpnTelCf4kT+eXxZ5cfp+15701BeBl4pKK5v8xi29D1k5ypybemayK09qQct7cfddFmGU/Evkxv62tnQ13tVH5vh8PMtYmebiuCXn+Zi4XbsIRknTk5PV3QT2/BznJ05dp3tdZIl8eW96SeTH6lqQtl1/J+ZXw6k5ZcU0iHE7WjN5Qa14oIbVjGLSg95q17rJckbmGWukueP0uKmoZRef3vIz3LgSTrKbV+yudr2+BvbN9GIYiO9TxcHzXq3dd/aO2ON045ZzbAg88hZPe1OtXoJU/6mH/jf9wR9BJp//wBMP/G/McavKOUqUEorm9CodvV9CJuLf2iGSf8ADfmcdCloyz/Wb7voxwtG41Etad0kQlZPuN3hLgApoYOTqQVJfeu+6xTow3muXGxs0qtOllTpytxvJNt+BnKt447XMLh1BWRepZGZDaS+4/FFlbRS9x+JxsrvGmtBskZq21Ffw5eKHf6vH7kvFGdVpclEilBcitLa6XuPxQxbVi/cfihqtbi16pEtOmUVteK9x+KB7aj/AC5eKGqbjSURTInt6K/hS8UR/tHD+VPxQ41OeMbSVxWjFj6SwX8Kf4kEvSaH8qf4kXhkzfJj9askYm0s6tTo7/GyY6XpHD+VLxRmYjaO/wAGm23LPW+nyOuGNnbj5MpZ6Wd3JPgxtiCG0Eobri2+dxv21fdfidHJZAq/bV91h9tX3WSxViayfcz16l7Ef6V9Dxl4xfdfidrT/wDkOiopfZ6mSS9uIkK0/Sqyp2S7dWUKafS98iX0dwcadOtSivZq3fNq0ZL5fQ5zF+m1GrWpTdCe5TbklvRu5WsixT9PqEKkpLDVO0lftR1V8/D6DXs360penTf2tPgoWXwf+TlKzzZr+km3oY2rGcKcoWjZptO+bZiSlcuk2WLJ1G0WysmTSrJ5WZREIx9g9U3oQRgTrCvmh/2J/eXgUVUdp6PYv1tJJ+1HJ9eTOSdC3FF/ZeM+zz3tU1mr2Oec5RvC6ruFElic9H0nh/Kl+JC/tZTX8GX4kcOGTvyn10iQyrKyOf8A2vp/yZ/iRBi/ShTg1ClKLeSbaaReGRzxc7UqJ8bJcEQN30GxfEdKWZ6XkKoXDuH03bUKcc7ICOWo2WuRLUi9WR7oCNcRYai72VhVTk3oUJMbYc4/EboQFshGDEKFuWsFjZ0ZqcJOMly49GVBUB6x6PbZhi6Sek17UeT8jWcTyHZG0p4aqpwemq5rkeq7L2hDE0ozi9V+kVE1RdmXc/oeOeuVl3Hs1ZdiX9L+h4nBLIzWpT9Ryo3HU4cixQpu/Qzt0xw2oTjZixL+IoJ95RUcy72mWHGn07xd0aNCrfoypCGXMk3d1mK7Y46WqryCnX4SGRvJW/5HrC20d/1yI3oSpKzaZGp7pL6qXukUlnnl3ENH7ylwsSwpx+613hQjd6d46vWjayzIukbSbdkNcOYirr/gTfT1Aikk9Fcr1KRowUZLLJ9SOrRazt+ZZWLjtlSiMZeqQK1SmdJXDLHSEQWwhtzAAAAAAACsQUBAAQBRAAAHIaKA9Ncrk1OXFldC7xFaManIVUk2VKDd8jRoR559wFepR5DN1Lq+ho1YPR8c7a5deBXnRSI0qTty/MhbJaj5ELCE3iei/wBWK1yeg1x07hUQqIqpsjU2SQxMoppWz6FQKLJ6asrc+JB9pl0+Yv2qXT55hTpvO471LtfUhlXb4L5jlipJWsrBEvqk1krfUjU/kNWJlZrLMZv9F8wJXHK2oyUXpyCNdpWsvmNlUbdwBCyQ1TFdV55IoaAXEAcmb3oztuWGqpSf7qT7S5f7jAuKqjQHtkaqnSck7pxf0PFoGxs70sxOHp+riqco5231J26KzRib7FIkU7MsQxLKW8KptGbG8c7GzCrFRb1ZSauytGu1yBYh9CcW75JWlTjZXJI0OLdrmfHaE1wj4PzEeOm9Uvn5k410/ri38PKGija3Es0N3evJZ9clY5uO1Jq3Zh4PzFnteo+EfB+Y41qebCN3G4iKyjZvojPV276vgZv2+fKPz8xy2lNcI+D8xxrP9ca0KlVpWRBJ5WKksfN6qPg/Ma8ZLlH5k41P64rJPSRnfa5cl8ySG0Zx4R+KfmONSeTFoXt0JoVr5PMyXtCb4R8H5jft0uUfn5k4Vf64tudFcErFPFYa+eRXpbXqRVkoNdU/MZLak37sPB+ZZjlGbnjUFSm0R2JqmMlLVR+CZA5HSON0AC4lyoUBLgAoCXFuAAJcAAAABUOihtw3gFYsY3G3F32FW8JTcnZLLmalF2XZs2vC5jUMZKnfdSz7x8dozSyUfB+ZNErcVPjxXy82VMRlf9eBUW2aqVlGFu5+ZG9pzfuw8H5mbK1LEvqbK7KtS7fIWpjpy1S+fmQyqN8iyJaUtUN2Obz6Ip7/AHD44hq2SLURAAFQAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAH//Z\n",
"text/html": [
"\n",
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.display import YouTubeVideo\n",
"YouTubeVideo('XOxxPcy5Gr4')\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Algorithms for sampling simple probability distributions\n",
"\n",
"Let's say I have the ability to sample a fair coin toss $X \\in \\{0,1\\}$ as many times as I want (equal probability of 0 and 1).\n",
"\n",
"#### Problem\n",
"\n",
"1. Using my access to coin tosses, how can I sample a real-valued random variable $U$ that is uniformly distributed on $[0,1]$? Assume I have to output this to $k$ values of precision.\n",
"1. Using coin tosses, how can I sample a discrete random variable $Z$ from $[n]$ that has distribution $P(Z = i) = p_i$, $i=1, \\ldots, n$? "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## More sampling\n",
"\n",
"Let's try a bigger class of distributions.\n",
"\n",
"#### Problem\n",
"\n",
"1. How can I sample a random variable $Z$ from the standard *exponential distribution*? A standard exponential has PDF $p(Z = x) = \\exp(-x)$ and CDF $P(Z \\leq x) = 1 - \\exp(-x)$. (*Hint*: can you somehow transform a random variable $U$ that is uniform on $[0,1]$ into an exponential distribution?)\n",
"1. Take any distribution with CDF $f(x) = P(Z \\leq x)$. Assume we can invert $f$, so for any $p \\in (0,1)$ we can find $f^{-1}(p)$. How can we sample from this distribution using only coin tosses?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Sampling can get very tricky\n",
"\n",
"Here are some canonical examples of tricky (and sometimes #P-hard!) sampling problems:\n",
"\n",
"1. Sampling a random spanning tree in a graph (this is possible!)\n",
"1. Sampling a random cycle in a graph (not easy)\n",
"1. Sampling a random Euler cyle in a graph (not easy)\n",
"1. Sampling a random perfect matching in a graph (we have good approximate method here!)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## The permanent of a matrix, and perfect matchings\n",
"\n",
"Assume we have an $n \\times n$ binary matrix $M$. Recall that the *permanent* of a matrix is defined as\n",
"$$\\text{perm}(M) = \\sum_{\\sigma \\in S_n} \\prod_{i=1}^n M_{i\\sigma(i)}$$\n",
"where $S_n$ is the set of all $n!$ permutations. \n",
"\n",
"**Note**: Computing the determinant of a matrix is easy! But computing the permanent is #P-hard!!!\n",
"\n",
"#### Problem\n",
"\n",
"Convince yourself that computing the permanent of binary matrix $M$ is equivalent to counting the number of perfect matchings in an undirected bipartite graph with $n$ nodes in each partition."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Sampling perfect matchings ==> Estimating the Permanent\n",
"\n",
"$$\\text{perm}(M) = \\sum_{\\sigma \\in S_n} \\prod_{i=1}^n M_{i\\sigma(i)}$$\n",
"\n",
"\n",
"#### Problem\n",
"\n",
"Assume you had some algorithm that can sample perfect matchings uniformly at random. How can you use that to estimate the permanent of a matrix?\n",
"\n",
"*Hint*: The trick is to construct a sequence of graphs, where each successive graph has one edge deleted. If $M'$ has one fewer edge than $M$, can you think of a way to estimate $\\frac{\\text{perm}(M)}{\\text{perm}(M')}$?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Markov Chains\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Markov Models\n",
"\n",
"A **Markov chain** is a series of random variables $X_1, \\dots, X_N$ satisfying the *Markov property*:\n",
"> The future is independent of the past, given the present.\n",
" $$\n",
" p(x_n | x_1, \\dots, x_{n-1}) = p(x_n | x_{n-1})\n",
" $$\n",
" \n",
"A chain is **stationary** if the transition probabilities do not change with time."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Markov Models: The \"graphical model\"\n",
"\n",
"We can represent a Markov chain using a \"directed graphical model\". The graph specifies how the sequence of random variables relate to each other, and in particular the independence assumptions.\n",
"\n",
"> This idea is complex and I'd recommend a course on graphical models to learn more"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true,
"jupyter": {
"outputs_hidden": true,
"source_hidden": true
},
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"@pgm_render\n",
"def pgm_markov_chain():\n",
" pgm = daft.PGM([6, 6], origin=[0, 0])\n",
"\n",
" # Nodes\n",
" pgm.add_node(daft.Node(\"x1\", r\"$X_n$\", 2, 2.5))\n",
" pgm.add_node(daft.Node(\"x2\", r\"$X_2$\", 3, 2.5))\n",
" pgm.add_node(daft.Node(\"ellipsis\", r\" . . . \", 3.7, 2.5, \n",
" offset=(0, 0), plot_params={\"ec\" : \"none\"}))\n",
" pgm.add_node(daft.Node(\"ellipsis_end\", r\"\", 3.7, 2.5, \n",
" offset=(0, 0), plot_params={\"ec\" : \"none\"}))\n",
" pgm.add_node(daft.Node(\"xN\", r\"$X_N$\", 4.5, 2.5))\n",
"\n",
" # Add in the edges.\n",
" pgm.add_edge(\"x1\", \"x2\", head_length=0.08)\n",
" pgm.add_edge(\"x2\", \"ellipsis\", head_length=0.08)\n",
" pgm.add_edge(\"ellipsis_end\", \"xN\", head_length=0.08)\n",
" \n",
" return pgm;"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"jupyter": {
"source_hidden": true
},
"slideshow": {
"slide_type": "-"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAB7CAYAAACGqB26AABA6ElEQVR4nO3dd1gU1/c/8PcWQDoi\nRAE1Knaxa7CXaOwQK9EYjRorURNjihpjiS35JJYYgy222DU2lKjRiC12UWIPtqgUUVFApGw5vz/4\nsj9nZ5e6u7ML5/U8+/DsZWfm7Ny9d8/euTMjIyJCMZKSkoLY2FjEx8cjPj4ecXFxiI+PR2JiIlQq\nFVQqFTQaDZRKJZRKJRwcHFC2bFn4+PjA19cXPj4+8PHxgZ+fH5ydnaV+O4wVaxkZGbr2mtNW4+Li\n8PjxY2RkZEClUkGtVkOhUECpVMLOzg5eXl6Cturr6ws/Pz94eHhI/XYYY8zmaLVaXc6k3xenpqZC\nrVZDpVJBJpPpcic3NzdR3uTr64ty5cpBLpdL/ZZMRil1AEWRlJSEixcvCh737t0zybplMhmqV6+O\nxo0b6x4NGzaEm5ubSdbPWEnz6tUrREdH48KFC7r2ev36dWi1WpOsv3z58oL22rhxY5QtW9Yk62aM\nseJAo9Hg33//FeRNly5dwsuXL02yfldXVzRs2FDQD1evXt1mE2eZLY0kq1QqnDhxAnv37kVERARi\nYmIsun2ZTIaAgAD06NEDwcHBeOutt2y24hkzNyJCdHQ0wsPDsXfvXkRFRZksIc6vihUrolu3bggO\nDkb79u1RqlQpi26fMcakFh8fj3379mHv3r2IjIw0WUKcX66urmjfvj2Cg4PRo0cPmxq8sPokOSUl\nBX/88QfCw8Pxxx9/IDk5WeqQdMqWLatLmDt16sRfwKzEU6lUiIyMRHh4OMLDw/Hw4UOpQ9JxdnZG\np06dEBwcjKCgIJQpU0bqkBhjzCyuXbuG3bt3Izw8HOfOnZM6HB2ZTIbAwEAEBQWhZ8+eqF27ttQh\n5cpqk+R//vkHS5cuxfr165GWllbg5T08PATzZcqVKwdHR0colUooFApoNBqo1WqkpaWJ5uGkpqYW\neHtlypTBsGHDMHr0aFSpUqXAyzNmy2JjY7Fy5UqsWLEC8fHxBV7eyclJ11Zz/rq4uOjmv2m1WqjV\namRkZODx48e69hoXF4ekpKQCb8/BwQEhISEIDQ1FYGAgZDJZgdfBGGPWJCMjA9u3b0dYWBjOnDlT\n4OUVCgXKlSsn6Is9PT1hZ2cHpTJ7dm7O/OSkpCRdHxwfH4+EhARoNJoCb7NFixYIDQ1F37594eDg\nUODlzc2qkuTMzEzs3LkTYWFhOHnyZL6WcXBwQP369UXzX5ycnAodR2pqKm7cuCGYs3P16lWo1eo8\nl5XJZOjSpQtCQ0PRtWtXKBSKQsfBmDUjIkRGRiIsLAy7d+/OVwcpl8tRq1YtQXsNCAiAm5tboRPV\njIwM3LlzRzTH7tWrV/lavmHDhggNDcWAAQP4ZF3GmM25d+8eli1bhlWrVuHZs2f5WsbX11d0zpWP\nj0+hp5BqNBrEx8fj0qVLgr44v4Mm3t7e+OijjzBq1ChUqlSpUDGYBVmBzMxM+vnnn6ls2bIEIM9H\n48aNaebMmXT+/HnKysqySIzp6el04sQJmjRpEtWuXTtfcVapUoXWr19PGo3GIjEyZglarZZ2795N\nAQEB+WoHlStXpk8//ZQOHz5ML1++tEiMarWa/vnnH/r++++pZcuWJJPJ8ozTw8ODZs+ebbEYGWOs\nKG7dukUhISH56t9cXV0pJCSENmzYQHFxcRaLMTY2ln777Tfq168fubi45BmnTCaj/v37U0xMjMVi\nzI2kSbJGo6ENGzZQ5cqVc91pCoWCunbtSsuWLaNHjx5JGbJOTEwMLViwgNq0aZNnpderV4/27dtH\nWq1W6rAZK5Jjx45R8+bN8/zMN23alObOnUtXr161is99YmIirVmzhnr16kX29va5xl6uXDkKCwuz\n2A9wxhgriNjYWBo5ciQpFIpc+zJfX18aN24cHTp0iDIzM6UOmzIyMujgwYP08ccfU7ly5XKNXalU\n0pgxYyya0BsiSZKs1WopIiKC6tWrl+eX1bRp0+jhw4dShJlvt27dogkTJpCHh0eu76d169b0999/\nSx0uYwUWHR1N3bp1y/Xz7eTkRCNHjqRLly5JHW6uEhMT6bvvvqNKlSrl+n78/f1p8+bNfCSIMWYV\nkpKS6KuvviJHR8dc+64OHTrQjh07rPqHflZWFm3fvp3at2+f5/fKlClT6MWLF5LEafEkOT4+nnr2\n7JnrTmnTpg1t3brVqivYkLS0NFq1ahU1atQo1/c3YsQISk5OljpcxvKUnp5OX375JcnlcqOf5xo1\natDixYsl68QKS61W0759+/JM/tu1a0d3796VOlzGWAm2c+dOeuONN4z2U25ubvTJJ5/QjRs3pA61\nwK5du0bjxo0jV1fXXAdNw8PDLR6bxZJkrVZLGzduJE9PT6M7oUWLFnT8+HFLhWQ2Wq2WwsPDc52z\nWaFCBTp48KDUoTJm1JkzZ6hmzZpGP8NVqlShjRs3FouR1qioKOrSpYvR9+rs7ExLliwpFu+VMWY7\nnj59SgMGDDDaNzk6OtKkSZMoKSlJ6lCL7NmzZ/TFF19QqVKljL7fDz74gJ49e2axmCySJOc1ehwQ\nEEDh4eFWMXfRlNRqNf3222/05ptv8qgysxl5jR6/8cYbtGTJEquY42ZqkZGRFBgYyKPKjDHJ5TZ6\nrFAoaNSoURQbGyt1mCb38OFDGj58uNHvIEuOKps9Sd6/f7/R0WNvb29au3YtqdVqc4chqYyMDFq4\ncKHRMzsrVqxI58+flzpMxujWrVtGr95ib29P06dPp9TUVKnDNCutVks7d+6kihUrGh1V3rhxo9Rh\nMsaKqfT0dBoyZIjRH+s9evSgW7duSR2m2d24cYO6du1qdD8MHz6cMjIyzBqD2ZJkrVZLP/74o9Ff\nAiEhIZSYmGiuzVule/fu0dtvv21wf5QqVYq/eJmkDhw4QO7u7gY/n02aNKGrV69KHaJFpaSk0KhR\no4x20JMmTSr2P/AZY5YVGxtLb731lsE+p3Tp0rR+/fpid9Q9N1qtltauXWv0u6l58+YUHx9vtu2b\nJUlOT0+nwYMHGx093r59uzk2axO0Wi0tXbrU6Kgyf/EyS9NqtTR//nyDP2jt7e1p7ty5pFKppA5T\nMocOHTI6qty9e3eeLsUYM4mzZ8+Sj4+Pwb4mKChI8suhSenRo0dGT7IuX748XbhwwSzbNXmSHBcX\nZ3ROX69evUrc6LEx9+7dM3qNZf7iZZaSnp5OH374ocHPYf369Uvc6LExKSkp9NFHHxncT7Vq1bKa\nC98zxmzT+vXrycHBQdS/uLi4lLjRY2NyRpWdnZ0NHo3fvHmzybdp0iQ5JiaGKlSoYPCLZPbs2VzJ\nerKysig0NNTg/qpXrx49fvxY6hBZMZaSkkJt27Y1Oh0qLS1N6hCtztKlS0mpVIr2l6enJ128eFHq\n8BhjNmjevHkG++EqVarwQIUB0dHRRi+I8MMPP5h0WyZLkm/cuGHwMIGzszPt2rXLVJsplox98daq\nVatEH15h5vPixQujd87jH7S5O3r0KJUpU0a039zd3en06dNSh8cYsxFarZamTZtmsB9u3749PX36\nVOoQrVZiYqLRo/GzZs0y2XZMkiTfuHHD4GVKKleuTP/8848pNlHsHT16lLy8vET7sFq1apwoM5N6\n8eIFNW3alH/QFsHdu3epbt26Bg+Nnjp1SurwGGM2YOrUqQaTvLFjx9rczdSkkJmZSaNHjza4D2fM\nmGGSbRQ5SY6JiTE4gtywYUOef1xAd+7cMXir3Nq1a/O+ZCaRmppKLVq0EH3GvLy8rP520tYmNTXV\n4NVq3N3dzXYSCWOseJg1a5bB5G7hwoVSh2ZzfvjhB4P78rvvvivyuouUJMfFxRmcgxwYGEjPnz8v\ncnAl0YMHD6hatWqifdqgQYNif31aZl5ZWVkGkzofHx+6du2a1OHZpFevXhm8jqenpyfdvHlT6vAY\nY1bop59+MpjULVu2TOrQbNaSJUsM7tOwsLAirbfQSXJ6errBq1gEBgbylRmKKC4ujqpWrSrat716\n9eLb4rJCGzNmjOgzVbZs2RJxUXpzysjIMHhpourVq/NgAWNM4ODBgwYvt7ly5UqpQ7N5S5cuFe1X\nhUJBf/31V6HXWagkWavV0qBBg0TBNGrUiL8UTOTBgwcGp15MmzZN6tCYDQoLCzM4xYLPnDaN9PR0\n6tixo2gfd+7cma97zhgjouw7mnp4eIj6iSVLlkgdWrGxcOFCg0f2bt++Xaj1FSpJ/vHHHw2epMfz\nZk0rJiaGSpcuLdrX27Ztkzo0ZkOOHj0qunpKqVKl+FboJvby5Utq1KiRqL1OnDhR6tAYYxJ78eIF\n1ahRw2wnmLH/7+uvvxbt5zp16lBKSkqB1yUjIkIBHDhwAN27d4dWq9WVubi44PTp0wgICCjIqlg+\n/PXXX+jcuTM0Go2uzNHREX///TcaNmwoYWTMFty7dw9NmzbFs2fPBOWbNm3CgAEDJIqq+Hr48CGa\nNGmCxMREQfm6deswePBgiaJijElJo9EgODgYf/zxh6C8T58+2LZtG+RyuUSRFU9arRa9e/fGnj17\nBOXBwcHYtWtXgfZ3gWomISEBAwcOFCTIALBhwwZOkM2kQ4cOWLhwoaAsPT0d/fr1Q1pamkRRMVug\nVqvRv39/UYI8efJkTpDNpEKFCti5cyfs7OwE5SNHjsTNmzcliooxJqUffvhBlCDXr18f69at4wTZ\nDORyOdavX486deoIysPDw0X5VJ7ryu8LiQijR49GUlKSoHz27Nl49913C7RRVjBjx47F8OHDBWV3\n7tzBlClTJIqI2YL58+fj3LlzgrKgoCDMnj1boohKhpYtW2LZsmWCsszMTAwbNkxwRIgxVvxdu3YN\n06dPF5R5e3tjz549cHZ2liiq4s/V1RXh4eHw9PQUlE+dOhW3bt3K93ryPd1i06ZNGDhwoKCsd+/e\n+P333yGTyfK9QVY4WVlZaNOmDc6ePSsoP3bsGNq0aSNRVMxaXb9+HQ0bNkRWVpaurFq1arhw4QLc\n3NwkjKzk+PjjjxEWFiYo+/HHHzFx4kSJImKMWZJarUaLFi1w/vx5XZlcLseRI0fQtm1bCSMrOf76\n6y+88847eD3Vbd68OU6cOAGFQpHn8vkaSU5ISMC4ceMEZd7e3li+fDknyBZib2+PdevWoVSpUoLy\nYcOG8bQLJqBWqzFkyBBBgiyTybBu3TpOkC3ohx9+QNWqVQVlBR3FYIzZrh9//FGQIAPAZ599xgmy\nBXXo0AHjx48XlJ0+fRqLFi3K1/J5JsnGplksXboUXl5e+Y+UFVmNGjVEh8p52gXTN3/+fIMdc/Pm\nzSWKqGRycnLCmjVrBAMJGRkZGDp0KE+7YKyYMzTNokaNGvj2228liqjkmjNnDvz9/QVl+R2wyHO6\nRUREBHr06CEoe++997Bly5ZChMqKSqPRoE2bNjh16pSuTCaT4dKlS6hfv76EkTFr8OjRI1SrVg0Z\nGRm6sho1auDSpUtwdHSUMLKSa8KECaJRi5UrV4rOM2CMFQ9EhA4dOiAyMlJXJpfLcfLkSR6skMiJ\nEyfQtm1bwbSLTp064eDBg7kul+tIslarxeTJkwVl3t7eWLJkSRFCZUWhUCiwevVqwbQLIsLXX38t\nYVTMWsycOVOQIMtkMqxZs4YTZAnNmTNHNO1ixowZSE9Plygixpg5HTp0SJAgA3w0T2qtW7cWTbv4\n888/ceTIkVyXyzVJ3rRpE65cuSIo+/HHH3mahcRq1KghmmIRERGBEydOSBQRswY3b97E6tWrBWWj\nRo3ijlliTk5OWLp0qaAsNjaWBxsYK4a0Wi0mTZokKCtfvjxPs7ACc+bMgY+Pj6Bs8uTJyG1ChdHp\nFpmZmahZsybu37+vKwsICMDly5fzdUYgM6+XL1/C399fcNOC5s2b4++//+aTKUuovn37YseOHbrn\nTk5OuH37tqhTYNJ45513cPjwYd3z0qVL4+7du/Dw8JAuKMaYSW3ZskV0HfpVq1Zh2LBhEkXEXrd8\n+XKMHj1aUPb777+jT58+Bl9vdCR5xYoVggQZAObNm8cJspVwcXHBN998Iyg7ffo09u7dK1FETErn\nz58XJMgA8Omnn3KCbEXmzZsneP78+XP873//kygaxpipqVQqTJ06VVBWs2ZNvtumFRk2bBiqVasm\nKPv666+hVqsNvt7gSHJmZibefPNNPH78WFfWsmVLnDhxgkcprUhWVhZq1qyJe/fu6crq16+PS5cu\ncT2VMN26dcP+/ft1zz09PXH37l24u7tLGBXTFxISgu3bt+ueOzo64uHDhyhTpoyEUTHGTGH16tX4\n6KOPBGU7d+5Er169JIqIGbJt2za89957grLffvsNgwYNEr3W4Ejyjh07BAkyAHz33XeceFkZe3t7\nzJo1S1AWHR0tuPIFK/7u3LkjSJCB7HlWnCBbn9mzZwuOxqWnp2Pt2rXSBcQYMwkiEp1nEBgYiJ49\ne0oTEDOqb9++aNSokaDsl19+Mfhag0my/l2i2rZti1atWpkoPGZKAwYMEF3/T7/+WPGmf1JY6dKl\nERoaKlE0LDfVq1cXjWAsXboUWq1WoogYY6Zw7tw5XLp0SVD29ddf8+CiFZLL5aJpMWfPnsXFixfF\nr9UviI6Oxt9//y0o+/jjj00cIjMVuVwumoS+fft2wQl9rPhKT08XXdFi6NChcHJykigilhf9HzB3\n7tzBoUOHJIqGMWYK+oNTb775Jrp16yZRNCwvQUFB8PPzE5TpDzgBBpJk/Rf5+Pjw4QIrN3ToUMF1\nk1UqFVatWiVhRMxStm7diufPnwvK9H80MevSokUL1KtXT1DGR38Ys11Pnz7F1q1bBWWjR4/mCx1Y\nMaVSiVGjRgnKNm3aJPo+FSTJKSkp2LBhg+AFI0eOhJ2dnZnCZKZQpkwZ9O/fX1C2bNkyvvVtCaCf\nXHXq1El05i6zLjKZTDSavG/fPvz3338SRcQYK4o1a9YgMzNT99ze3l50Ah+zPiNGjIBSqdQ9T09P\nx7p16wSvESTJERERSEtL0z1XKBQYMWKEmcNkpqD/pfvgwQOcPXtWomiYJdy9exfnz58XlPFcZNsw\ncOBAuLq66p5rtVrRJfwYY7ZBfxQ5JCQE3t7eEkXD8qtcuXKi6yPr16UgSQ4PDxf8s3PnzqI5G8w6\nNW3aFHXq1BGU6dcnK170r4nt5eWF7t27SxQNKwgXFxeEhIQIyri9MmZ7Hj16JDrha+jQoRJFwwpK\nv67Onj2LhIQE3XNdkpyVlSW6jNS7775r5vAKbu7cuZDL5fl6NGnSxOA6QkND81xWoVDg999/t/C7\nKxr9+uIv3eJNv3579OghOHRkTbjdium315MnT+LZs2cSRcMYK4x9+/YJnnt4eKB169YSRZON+9v8\na9euneCoHhEhIiICrxcQEdHhw4cJgOARGxtL1iY9PZ3u3btHFy5coLCwMCpTpgzJ5XLdw8HBgaZN\nm0YXLlyg1NRUg+tITU2lI0eOkL+/v2DZunXr0pw5c+jo0aMUExNj4XdWdGfOnBHV4b///it1WMwM\nnj9/TkqlUlDXO3fulDoso7jdiqWlpZGjo6OgDtevXy91WIyxAujataugDb///vtSh8T9bQH169dP\nUIfBwcG6/+mS5PHjxwte1LRpU0mCLaiVK1eSTCbTVWCPHj3yvey4ceNIJpORp6cnrVq1yoxRWoZG\no6GyZcsK6nH+/PlSh8XMYPPmzYJ6dnBwMNrZWSNut9mCg4MF9divXz+pQ2KM5VNqairZ29sL2vCW\nLVukDkuE+9vcrV+/XlCHjo6OlJaWRkSvJcn+/v6CF3377beSBVwQKSkp5OzsTHK5nGQyGTk6OtKL\nFy/yXC49PZ0qVqxIfn5+dO3aNQtEahkfffSRoB7ffvttqUNiZvDBBx8I6rlbt25Sh1Qg3G6zrVy5\nUlCPrq6upFKppA6LMZYPu3fvFrRfpVKZr37M0ri/zd3Tp09JLpcL6jIiIoKIiOQA8OzZM9y5c0cw\nT8NWTgBydXVFnz59QEQAgMzMTGzcuDHP5QYOHIiMjAwcO3YMtWvXNneYFqNfb+fPn+e7eRVD+lcu\nsZX2moPbbTb9ektNTcXNmzclioYxVhD6/XCrVq3g7u4uUTTGcX+buzJlyqB58+aCspy6lQMQnZlZ\nqlQp0cXurVnO2YkymQxEhDVr1uT6+i+++AL79+/H7t27Rbd0tnWBgYGC56mpqYiJiZEoGmYOycnJ\nojrVr3dbwO02+2ZNFStWFJQZujUqY8z66LdVa+6Hub/NnX7d5dStHAAuXLgg+GeDBg2s9ix5Q9q1\na4fKlSvrnkdFReHq1asGX7tixQosWLAAq1evFv1yKA58fX3h4+MjKOMv3eIlKipK8Nze3h4BAQES\nRVN43G6zNW7cWPBcvz9mjFkfIhK1VWNXirAG3N/mTr/uBEmyfhKl32nbgiFDhoCIIJPJAACrV68W\nvebPP//E2LFjMWPGDNEd6ooT/frjJLl40a/PunXrwsHBQaJoiobbLbdXxmzRf//9h6SkJEGZtedO\n3N8ap193CQkJiIuLK15JslyefdlnIsLGjRuhVqt1/79y5QpCQkLQv39/fPPNN1KFaRH8pWubXr16\nla/XFYf2moPbrXgE4/Lly4J9IKX8fiYZKy4K2w+XLl0alSpVMkNEpsP9rXFVq1YVXC8ZyK5jeXJy\nMv777z/BP2zxS7dChQp4++23dRPTnz59qrsjWXx8PHr06IH69etj1apVUoZpEfr1Fx0dLVEkrCCa\nNWuGHTt26D7DxujXpy221xzcbsX1l56ejtu3b0sUjdDUqVPx+eef4+XLl1KHwphFfPXVV5g0aVKe\nybKhfjhnhNZacX9rnFwuR6NGjQRl0dHRkMfGxopeXL16dUvFZVKvT0wHsg8lpKenIygoCKVKlcKu\nXbtgZ2cnZYgWUaNGDcHzFy9eIC0tTaJoWH7Fxsaib9++CA4OFv1w1X/d6/Tr29aU9Hbr5eUFT09P\nQZmhflkKarUa8+fPR506dUR3FmOsOFKpVPj+++9Rp04d0V2IX2er/XBJ729zo1+HsbGxkMfFxQkK\nPT09UapUKUvGZTK9e/eGh4cHgOxDCQcPHsS7776L+/fvIyIiQvRFVFzpn7gHZP9KZLZh3759qF27\nNubPny867J6WloaUlBRBmaH6tiXcbsV1qN8vS+3BgwcICgpC3759rSaBZ8yc7t+/j27duqF///4G\nvz/126ivr6+lQisS7m+N06/DuLg4yPUr35a/cB0cHNC/f3/dxHS1Wo3jx49j586dqFq1qtThWYyr\nqytcXFwEZZwk25ZXr17h888/R9OmTXHu3DlduaF6tOU2C3C7BcSds7W21x07dqBWrVpYsmQJNBqN\n1OEwZnZbt25FrVq1sGzZMsE9B2w1d+L+1jhDgxWikWRb+TVkzLBhwwTPK1WqhDZt2kgUjXQM/SJi\ntufy5cto1qwZxo0bh+TkZFE9uri4iE42sEUlvd1a+0jy61JTUzFu3Dg0b94cly9fljocxswuOTkZ\nY8aMQcuWLXHlyhUAtjuSDHB/a4yhwYpiNZIMZJ8pXqdOHd3zmJgYnDlzRsKIpKFfj9Y6MsXyRkRY\nsmQJateujd27dwv+Z+vtNUdJb7e22F7Pnz+PJk2a8Il9rMQ4c+YMGjVqhC+++AJPnjwR/M+W+uKS\n3t8ao1+HCQkJkOtXdLly5SwZk8k9ePAAz549E5SVxDM19esxMTFRokiYqcTFxWHhwoWCMltvrzlK\neru11faq0Wj4xD5WoqjVavz444+iclvqi0t6f2uMfh1qNBrIs7KyBIWOjo6WjMmkXr58iR49euCN\nN96Avb297vaL27ZtQ3p6utThWZR+Pc6bNw8ymYwfVvzQvzB9ftjqSbav43Yrbq9Hjx6V/PMok8nw\n888/5yt+PrGveJP6c2jJx/Llywu1j2wld+L+1jhDdShXqVSCAlu6HfXrtFotQkJC8OLFC+zfvx9B\nQUG6awG+fPkS27ZtkzhCy7LVemQFc+bMGcGJfbbGnO1Wq9VixYoVaN26NRo2bIiKFSuiefPmWLRo\nkdVdErG4tFc+sY+VVGvWrBGc2GeNLJEnPXr0CK1atUKTJk3g5eUFuVwOuVyOL7/80ugyV65cQaNG\njVCqVCnd6+3t7REQEIBFixYVOab8MtgP9+jRgwDoHnPnziVb9PHHH5OLiwtdvnyZiIj2799PMpmM\n5HI5yeVyatOmjcQRWtaoUaME9cqP4vtwdnbWfe5tjbna7atXr6hjx4706aefUlpaGhERZWVl0Zw5\nc0gmk5G/vz9dunTJVG+jyNasWSP558jUj2+//Vbq3cpMROrPkq08vvvuO6mrKleWzpO++uor3bo9\nPT0pPT0919fHxcWRh4cHde7cmZ48eWLSWPLj5cuX4nrt2bOnoGDWrFkWD6yoFi1aREqlksLDw3Vl\nWq2WypcvT3K5XPchiImJkTBKyxo+fLjkHQY/zP/w8vKi27dvS/1xKxRzttvBgwdTaGiowf9NnjyZ\nZDIZ+fn5UXx8fKHjN6Vff/1V8s+SqR5+fn60c+dO0mq1Uu9WZiJSf6Zs4bFt2zapqylXUuRJzZs3\np/fee0+37lWrVuW5TMWKFenhw4cmi6EgkpOTxXXbt29fQcH06dMlCa6w9u7dSwqFghYsWCD639Sp\nU0kmk+kqaPLkyRJEKI0hQ4ZI3mnww/yPTp06Sf1RKxRzttt///2XHB0d6enTpwb/n5ycTA4ODiSX\ny+mTTz4pTPgmt2zZMsk/S0V9yGQyGj9+PCUnJ0u9O5mJyeVyyT9f1v5ISUmRupqMkiJPSk9PJycn\nJ0pMTCRvb2+Sy+XUpEmTXJd5+PAhBQQEmGT7hfHs2TNRvcodHBzwOmubq5eby5cv4/3338eIESMw\nYcIE0f+HDh2qm4xPRPjtt99082+KO/37zk+ZMgVExA8rfhTmTke21F5zmLvdHjlyBBkZGahZs6bB\nKy64ubkhMDAQRISNGzcW6b2Yin57bd++veSfRyLCuHHj8hV/w4YNcfbsWfz0009wc3Mzxy5iEtJo\nNJJ/Fi31GDVqVKH2kbX2xVLlSadPn0bt2rXh7e2NESNGgIgQFRWV6zk0x44dk/R6zfr9MADIy5Yt\nKyhISEiwVDxFEhcXh+DgYDRv3hy//PKLwddUqVIFbdq00VV4fHw8/vjjD0uGKRn966za0uVpmGF1\n6tTBrFmzBGW20l5zWKLdpqamAgCSkpKwcuVKg68pX7687jX6t/mWgq22V2dnZyxYsADnzp1D06ZN\npQ6HMbNzcXHBwoULIZfLBeXW2BdLmSedOHEC7dq1AwCMGTMGCoUCABAWFmZ0mePHj0uaJOv3w3Z2\ndpDb0p2ecqSnpyMoKAhubm7Yvn276MP6upw7y8hkMgDA6tWrLRKj1PTr0ZYudM6ESpUqhXnz5iEq\nKgrt27cX/C8uLs4kv/otwVLtdsCAAahbty68vLwwcuRIo7HksIb9Z4vtNSgoCNevX8eECROKzdU5\nGMtN7969cePGDXz66ad44403BP+zttxJ6jzp+PHjaNu2LYDsQYng4GAQEbZv347nz58bXObEiRO6\nZaSgX4flypUTJ8nWfqcnIsKAAQPw6NEjRERE5Hlor2/fvrrXEBEiIiLw9OnTfG9v0aJFaNOmDWrV\nqiW4kcOpU6cQEhKCDh06oE6dOmjVqhVOnjxZuDdlYkRks/eVZ0KdOnXC1atXMWnSJNjb24vqMT09\n3SpGQvNiyXbr5+eH6OhoPH78GN27dzf4mn/++Uf3Wnd39wK8E/Owpfbq5+eHnTt3Ys+ePahYsaLU\n4TBmdhUqVMCePXuwY8cO3VEoa86dLJ0n6VOr1Th37pxgVDhn6lZGRobBG5c8efIEGo1G0qNo+nXo\n6+sLuf69qq3t15C+CRMm4M8//8SePXvw5ptv5vl6R0dHvPfee7rRIpVKhfXr1+drW5GRkbh+/TqO\nHz+OYcOGYeLEidi8eTPmzp2LLVu2YNmyZfjrr79w5coVqFQqdOvWzSoupJ+SkiKaW2NL95VnQNmy\nZbF582YcOHAA/v7+unJDyZM1dc7GWLLd5iU6Ohp3796FTCbD4MGDTbLOojLUOVsbmUyG8ePH4/r1\n6+jVq5du1Imx4koul+Ozzz7D9evXERwcLPifNedOUve3Fy9eRLVq1QTJebt27VCnTh0QkcEbtpw4\ncULSqRaAuA59fX3FI8nJyckGJy9bg48//hiLFy/GDz/8gGbNmuV7udcPJRBRvg8lzJ8/H5999hkA\n6A5VfPnllyAiLF68WHeilVwux9tvv420tDTs2bOnIG/JLAwlTdY8MsWERo0ahRs3bqB///6iRMTR\n0REeHh6CMmvqnA2xdLvNy7x58wAA3t7e+Pzzz02yzqKy9ukW1n5i3rFjx3Q3IZDL5br5j+bQrl07\nwba+/fZbs2zHku+JiTVp0gTnz5/H/Pnz4eLiIvq/tU5VtYb+9vjx47r5yPqxAcDdu3dx8OBB0TJS\nJ8mGjujJcw4dvO7GjRuWiilf7t+/j27dumHp0qWQyWQYOHBggZYPDAwUHBa8fv16nncpe/XqFR49\neoSaNWsCyL4jDADUrl0bX3/9tcHXA9mHEqR2/fp1wXMvL69icfvi4i4gIAAnT57EsmXLULp0aaOv\nq1ChguC5tbXXHFK027wcOnQI27dvh1KpxMaNG0U/OKQQFxeH5ORkQZl+HUvFxcXFpk7Ms8Totv5t\njC2xPWY5rq6uWLx4Mc6cOYNGjRoZfZ219cPW1N++Ph/5dYMHD9ZNb9M/gc/YMpaknztVqFABICKq\nWrWq4Lpwy5cvJym9ePGCoqKiaOvWrTRgwABycnISXMfv888/z9cFrxMTE+nEiRM0cuRI3fI5D39/\nf1q+fDldvHiR7t27RyqVSrBsTEwMLVmyRPe8UqVKJJfL6fTp0wa39fbbb5NcLqedO3cW7c2bwJQp\nUwT12bFjR6lDYvmgVqvz9bpBgwYJ6nfo0KFmjix/rKHd5ub58+fk5+dHSqWSNmzYUJS3alLh4eGC\n+nR1dSWNRiN1WESU/8+k1I4ePar7nOX8NZd27doJPtczZ840y3Ys+Z7Y/5ffz/yePXskbbfW2t9q\ntVoqU6YMPX/+3OD/P/30U5LJZKRUKunBgwdElH3t+ipVqhTo/ZuaSqUiJycnQZ3u37+fQETUv39/\nwT9GjBghabA1a9bU3crQ2MPe3j7Xu7KoVCpycXHJcz05jx07dhhd1927d0kmk5G7u7vBRvDy5Usq\nVaoUKRQKozcwsKROnToJ6vOrr76SOiRmQj/99JOgfuvWrSt1SERkfe32dRqNhrp06UJ2dna0efNm\nU71lk5g+fbqgPk19a9iSICehfD1JMBdLJ8mWeE+s4GJjY0U3nrh586bFtm+t/W10dDQ1aNDA6P9j\nYmJ065syZQoREUVERNDgwYMLvhNM6MqVK6L6TExMJCUANG7cGFu2bNENMV+8eNFUo9eFYorDFkql\nUne91KKKjIwEALRu3drgZVQOHjyIzMxMtGjRAmXKlDHJNguLiET117hxY4miYeagX5/Xr19Heno6\nHB0dJYoom7W129eNGzcOx48fx65du4xe8UIq+u21SZMmEkVi216flmDuKQqWmgJhyffECsbX1xfl\nypUTXB/54sWLqFGjhkW2b639rbH5yDmqVq2Kzp0748CBA1i1ahVmzpxpFSft6ffDFSpUgLe3N+SA\n+Ev3ypUryMzMtFx0Vi4yMhIymUx0jdocO3bsgEwmQ79+/SwcmdiDBw/w7NkzQRknycVLgwYNBD/W\nNBoNoqOjJYzIui1YsAC//fYbDhw4IEqQ//vvP2i1Wokiy3bhwgXBc26vBde2bVtoNBrdQ61Wm21b\nkZGRgm1NmzbNLNux5HtihaPfVqUeYLQG+bnW8dixYwFkX/Zt27ZtVjEf2djgohyAaHK6SqVCVFSU\nhUKzfjkjyYaS5NTUVOzevRsKhQL9+/cHADx8+BCLFi2yZIg6p0+fFjwvXbo0KleuLEkszDycnZ11\nJ5Tm0K93lm3Xrl2YOXMmIiIi0Lp1a9H/e/bsCY1GI0Fk2e7fvy+6UxcnyYzZBv22yv0wcPLkyTwT\n3q5du+oubbpw4ULExcWhatWqlgjPKP26EyTJ7u7uqF27tuAF+/bts1Bo1u3ff/9FXFwcPDw80LBh\nQ9H/d+zYgfT0dHTq1Ak5t/jesmWLZDco0K+3wMBAPkxXDDVv3lzwnNur2Llz5zBy5Ejs2bPH4KG8\nrKwsqFQq2NnZSRBdNv168/T0RLVq1SSKhjFWEPr98JkzZ/DkyROJopHe7du34eXllevVmYDsqUNj\nxowBESEqKqpAl6ozh/j4eNERvZyYdMdse/ToIXhBeHi4BUKzfjmjyMbmy1y6dAkymQy9e/cGkH0J\nuK1bt0oy9UKtVovuuR4UFGTxOJj56bfXY8eOGb3VZ0l0//599OzZE9OnT4ePjw9u3boleqxfvx5V\nqlSRNE79frZ79+653j6WMWY92rVrBycnJ91z+r+71ZVUBw4cyPfc4mHDhun2ndTzkfUHK1xdXXUx\n6Xpj/bvJXL16Fffu3bNAeNYtr/nIAQEBALJvTEBE+Oyzz/DJJ58YvPi4uf3999+iRImT5OLpnXfe\ngYODg+65RqPBgQMHJIzIeiQnJ6Nbt25ISEjA+PHjUatWLYOPESNGiKatWFJKSgqOHj0qKNPvhxlj\n1qtUqVLo3LmzoGzv3r0SRSOt06dPY86cOXj06FG+Xu/h4aG7lrPU85H166xr166wt7cH8FqS3KxZ\nM3h5eeW6YEnk7OwMf39/hISEGPz/8OHDMX78eHzxxRdo1aoVateujUGDBlk4ymz6o1INGza0mpsS\nMNNydnZGx44dBWV89CdbWFgYbt26Jbrpg/5DLpdLmiQfPHgQKpVK99ze3l70hcsYs276P2wPHjxo\nFTcVs5Q+ffqgUqVKaNWqFRITExEeHg4fHx906dIlz2XHjRsHPz8/0XRfS3r16hUOHTokKHu9TmVE\n/3ezbgBDhw7F2rVrdf9s3749jhw5Yv4oWZEREapVq4Y7d+7oyqZPn44ZM2ZIFxQzq+XLl2P06NG6\n525ubkhMTBSMMDPrNXDgQGzatEn3POeySIwx2/HkyROULVsWr6VS2Ldvn9VdapIZtnv3bvTq1Uv3\nXKFQIDExEZ6engBeG0kGxL+IIiMjBUkXs17Hjh0T1RUfui3e9Oclp6SkYMeOHRJFwwoiKSkJO3fu\nFJRxe2XM9nh7e6NFixaCslWrVkkUDSso/bpq3bq1LkEG9JLkLl26iM5KXLZsmRnDY6aifx/0GjVq\nGLwaBys+/Pz8RBdt1/8cMOu0du1awSFZOzs79OnTR8KIGGOF9f777wue79mzJ99zc5l07t+/LzrR\ncsCAAYLngiTZ0dERQ4cOFbxg9erVSE9PN1OIzBTi4uKwa9cuQVloaChf+q0ECA0NFTz/+++/+cYi\nVk6r1WLp0qWCsr59++ouIckYsy0ffPCB4GR9rVaLFStWSBgRy4/ly5cLpsm4ubnpTibMIbrW0Otz\nHIHsw4Jbt241U4jMFFauXCm4G5OTkxMGDx4sYUTMUnr27Ily5coJyvQTMGZdDh8+jNu3bwvK9H/s\nMMZsh5ubm+iE/RUrViArK0uiiFheMjIy8OuvvwrKhgwZAmdnZ0GZKEmuVq0aOnXqJCj75ZdfBNk2\nsx4qlUr0i3XgwIHw8PCQJiBmUXZ2dhg5cqSgbMOGDUhOTpYoIpaXX375RfC8bt26aNmypUTRMMZM\nYcyYMYLnjx8/Fp13wKzH9u3b8fTpU0GZfh0CBpJkQDyqceHCBT7r2kqtWrUKcXFxgjIelSpZRo4c\nCYVCoXuelpaGBQsWSBgRM+bSpUuiS/Xx1CjGbF/dunXRunVrQdns2bMlve09M0ytVmPOnDmCsrff\nftvgJUENJsndu3dHpUqVBGWTJ0+GVqs1XZSsyNLS0jBz5kxBWcuWLdGgQQNpAmKS8PPz093xMcf8\n+fPx+PFjiSJixkyePFnw3N3dXTQHjjFmm8aOHSt4fu3aNWzYsEGiaJgxa9euxa1btwRl+nWXw2CS\nrFQqRdfXjY6OxpYtW0wTITOJn376CQkJCYKy2bNnSxQNk9KMGTMEtzNOS0vjz4KViYyMxMGDBwVl\nkyZNgqurq0QRMcZMqW/fvqhXr56gbNq0aSXq5iLWLj09XZTfNm7cGO+++67B1xtMkoHsszXr1Kkj\nKPvmm294IrqVePbsGb7//ntBWefOnUWXBGMlQ+3atfHhhx8KypYvX467d+9KFBF7HRFh0qRJgjIf\nHx+MHz9eoogYY6Yml8sxb948QdmDBw/4UrpWZMmSJYiNjRWUzZs3TzDI9DqjSbJCocDcuXMFZXfv\n3sXy5ctNECYrqnnz5iElJUVUxkquGTNmCO62p1KpMHXqVAkjYjl27tyJc+fOCcqmT58OJycniSJi\njJlD165d0aZNG0HZ7Nmz8eLFC2kCYjpJSUmiPKlDhw545513jC4juC21PiJCq1atcOrUKV2Zq6sr\nrl69iooVK5ogZFYYFy5cQLNmzQQnBAwYMEBwi1tWMk2cOFF00t4ff/yBrl27ShQRe/78OerUqYP4\n+HhdWbVq1XDt2jXY2dlJGBljzBxOnz4tugvfsGHD+E58Ehs8eDDWr18vKDt37hyaNm1qdJlck2QA\nOHHihOhXUadOnXDgwAE+I1sCmZmZaNKkCa5evaorUyqVuHHjBqpWrSphZMwaPH36FP7+/oKjDH5+\nfrh27Rrc3d0ljKzkGjJkCNatWyco27p1K0JCQiSKiDFmbr169cLu3bsFZfv370eXLl2kCaiE27t3\nL4KDgwVlffv2xfbt23Ndzuh0ixytW7cW3YXvzz//5F9EEpk1a5YgQQaAqVOncoLMAABeXl744Ycf\nBGWxsbH47LPPJIqoZIuIiBAlyF27dkW/fv0kiogxZgmLFi0SnZQ7YsQIvoa9BJ4/f45Ro0YJytzd\n3bFw4cI8l81zJBkAXrx4gYCAAMFkZ552YXmGplk0aNAA586d48O2TIeI0LlzZxw6dEhQztMuLMvQ\nNAs3Nzdcu3YN5cuXlzAyxpglrFixQpSc8bQLyzM0zWL16tWiAWBD8hxJBgAPDw+sXLlSUJaamooh\nQ4YIbofMzCc1NRUffvihIEFWKpVYs2YNJ8hMQCaT4ddffxWNYgwfPpyvnWwhRIQxY8YIEmQAWLhw\nISfIjJUQI0aMQMeOHQVlq1evxq5duySKqOTZvn27KEHu2rUrhgwZkq/l85Uk56xUP+uOjIzExIkT\n87sKVkharRaDBg3C9evXBeVTp07lG4cwgypWrIj58+cLyuLi4tCnTx9kZmZKFFXJ8f3332Pr1q2C\nMkN9KGOs+DI2YDF48GDRtElmepcvXxYlw+7u7lixYkX+z6mjAnj+/DmVL1+eAAgeK1euLMhqWAF9\n8803on3eoEEDysrKkjo0ZsW0Wi116tRJ9Nn56KOPSKvVSh1esbV3716SyWSCfe7u7k4PHz6UOjTG\nmARWrFgh6ocrV65MT58+lTq0Yuvx48dUsWJF0X5fvXp1gdaTrznJrzt79izatm0rGI2ys7PDkSNH\n0KpVq4KsiuXDtm3b8N577wnKPD09ce7cOfj7+0sUFbMViYmJaNq0KR48eCAoX7x4McaNGydRVMXX\n9evX0axZM6SmpurKZDIZ9u3bh27dukkYGWNMKkSEDz/8UHTYv3379jh48CBPmTSxrKwsdOzYESdO\nnBCUDxs2DL/++muBrsyW7+kWOQIDA7FixQpBmUqlQp8+fXD//v2Cro7lIioqSnSoQKFQYPv27Zwg\ns3x54403sGfPHtFNKyZMmCA6sY8VzbNnzxAcHCxIkAHgu+++4wSZsRJMJpNhxYoVeOuttwTlkZGR\n+OSTT1DAsUqWCyLC2LFjRQlyixYtEBYWVvBLFxd2KPvzzz8XDWP7+/vTo0ePCrtK9pqrV6+Sl5eX\naB///PPPUofGbNDvv/8u+iw5OzvTyZMnpQ6tWHj+/Dk1atRItI8HDhzIU1sYY0REFBsbS76+vqJ+\nYurUqVKHVixotVr68ssvRfu3fPnylJCQUKh1FjpJVqvV1KVLF1EwNWrUoPj4+MKulhHRzZs3qWzZ\nsqJ9O2LECP7CZYU2bdo00WfK1dWVzp49K3VoNi05OZkCAwNF+7ZJkyb06tUrqcNjjFmRc+fOkYOD\ng6i/mDVrltSh2bzp06eL9qujoyNdvHix0OssdJJMlD16UqdOHYOJMo8oF87Vq1cNJsjt27enzMxM\nqcNjNkyj0VBISIjos+Xm5sYjyoWUlJREb731lmifVqhQgftAxphBmzdvFp3cC4CmT5/OA2GFoNVq\nacqUKaL9KZfLadu2bUVad5GSZCKihIQEqlGjhii4KlWq0J07d4q6+hLl4sWLBqdYtGjRglJTU6UO\njxUDmZmZFBQUZHDqxeHDh6UOz6Y8fvyYGjRoINqXPj4+FBMTI3V4jDErtnLlSlHfAYC++OIL0mg0\nUodnMzQaDU2YMMHgvly7dm2R11/kJJmI6NGjR+Tv7y8KsEyZMhQZGWmKTRR727ZtI0dHR9E+bNq0\nKb148ULq8FgxkpGRYXCqlFKppLCwMKnDswlRUVFUoUIF0T5844036Pr161KHxxizAb/88ovB5K53\n7948MJYPKSkp9O677xrch8uXLzfJNkySJBNlJ8o1a9bkL94C0mg0Bq+DnDOCzAkyM4f09HSDI8oA\naPTo0Ty1JxfGftD6+vrSjRs3pA6PMWZDli1bZrAfrlevHt27d0/q8KzWnTt3DE73lclktGrVKpNt\nx2RJMlH21Iu6devyF28+paSkUM+ePQ3ur/bt2/MvSWZWmZmZ1K9fP4Ofv7Zt21JiYqLUIVoVjUZD\nU6dONbi/3nzzTbp9+7bUITLGbNDatWtJoVCI+hUvLy86evSo1OFZnb/++os8PT1F+0uhUNCGDRtM\nui2TJslE2Wd6GxuhatGiBc/V+z8XLlyg2rVrG9xPw4cP5x8UzCI0Gg19/fXXBj+HFStWpCNHjkgd\nolWIjY2lrl27GtxPrVq1osePH0sdImPMhh06dIhKly5t8Gj8vHnzSKVSSR2i5FQqFc2ePdvgDwpP\nT0+zfF+ZPEkmyv7iNXSmIZB9OY5FixaV2InpGRkZNHXqVIOVrFAo6Oeff+azW5nFbd68mUqVKmWw\nzYaGhpbYoxparZbWrVtHHh4e/IOWMWZWMTExVKtWLYN9zVtvvUXXrl2TOkTJXLlyhRo3bmxw3wQE\nBJjtQhFmSZJzbNq0yegXb+vWrUvcqPKFCxcoICDA4P7w9PSkv/76S+oQWQl24cIF8vPzM/j5rFy5\ncokbVY6NjaUePXoY3B/8g5YxZg7JycnUvXt3g/2Ovb09fffddyVqVDln9NjOzs7gPnn33XcpJSXF\nbNs3a5JMRHT+/HmqVKmS0VHlOXPm0MuXL80dhqSePXtGEydONDh6DIAaNGjA8xmZVYiPj6c2bdoY\n/JzmjJzGxsZKHaZZZWZm0s8//2x09Lhs2bJ8uTzGmNmo1Wr6+uuvDV5LGci+UdHx48elDtPsjh49\navBOpkD2NZCnTZtm9lkJZk+SiYhSU1MpNDTU6BdvuXLlKCwsjLKysiwRjsWkpaXR3Llzyd3d3eD7\nViqVNGPGDD5cy6yKRqOhRYsWGbyCQ86P20mTJtHz58+lDtWkNBoNbdiwgSpXrmy0r3r//ffp6dOn\nUofKGCsBTp06RdWrVzfaH3Xv3p3++ecfqcM0uUuXLhm8TGnOo1atWnTmzBmLxGKRJDnHkSNHjI4q\nA6CqVavSli1bbH6+clZWFi1dupR8fHyMvtf69evTpUuXpA6VMaNiYmKodevWRj/DpUuXpv/97382\nf+tlrVZLERERVK9ePaPvtWzZsrRr1y6pQ2WMlTCvXr2izz//3Oioskwmo0GDBtHdu3elDrXI7ty5\nQ++//77Rflgul9NXX31F6enpFovJokkyUfao8scff2x0JwCgatWq0YIFCygpKcnS4RVJXFwczZw5\n0+i8zpzR4+nTp/PoMbMJeY0qA9mXKfrqq69srpNOTU2l5cuXU/369XPtjwYMGMCjx4wxSZ06dcrg\n3Y1zHgqFgvr27UtHjhyxqXMltFotHT58mHr37m10SqqlR49fZ/EkOUdUVBR17tw51y8nR0dH+uij\nj+jixYtShZknrVZLR48epZCQEFIqlbm+n/79+5e4kxVZ8fDo0SMaMWJErp2YTCaj7t27U0REhFUf\nDbp+/TqNGzeO3Nzccm2vLVq0KBHz/hhjtiEzM5MWL15M3t7eufZdNWvWpMWLF1v1zcieP39OixYt\nyjXxzzmK98svv0g2sChZkpzjyJEj9NZbb+W6k4Dsu8988803dP78ecm/gFUqFR07dowmTpxIVatW\nzTP2zp07W3Wiz1h+3bx5k/r27ZvnZ75ChQo0duxYOnjwoORHTbRaLf3zzz80Z84cCgwMzDP2OnXq\n0J49e2xqNIYxVnKkpKTQzJkzydXVNde+zMnJiXr37k1r1661iptDPX78mFavXk09e/bM9egkAHJz\nc6PZs2dLfmEHGRERJEZE2L17N6ZNm4arV6/m+XofHx8EBQWhe/fuCAwMRNmyZc0e36NHj3D69Gns\n27cPERERSEpKynO5Fi1aYPbs2Wjfvr1Z42PM0s6fP4+pU6fizz//zPO1rq6u6NKlC4KCgtCyZUtU\nrlwZMpnMrPElJSXh3Llz2L9/P8LDw3H//v08l6lSpQq++eYbDBo0CAqFwqzxMcZYUT158gRz587F\nihUr8OrVq1xfK5PJ0KJFCwQHB6NDhw4ICAiAg4ODWePLzMzElStXcPjwYezduxenT59GXimns7Mz\nRo0ahSlTpqBMmTJmjS8/rCJJzkFEOHHiBMLCwrBjxw6o1ep8Lefn54fGjRvrHjVq1ICvry+cnJwK\nHENqairi4uJw/fp1XLx4Ufd48uRJvpZ3dHTEwIEDMWbMGDRq1KjA22fMlty4cQPLli3D2rVrkZKS\nkq9lSpcujUaNGunaa0BAAHx9feHu7l7g5DkjIwPx8fGIiYkRtNf8JMVA9hdH9+7dERoais6dO0Mu\nlxdo+4wxJrUXL17gt99+Q1hYGG7dupWvZezs7FC3bl1dP9ywYUNUrFgR3t7eBR4k0Gg0SExMxIMH\nD3Dp0iVdP3z16lWoVKp8raNWrVoIDQ3FoEGD4O7uXqDtm5NVJcmvS0hIwK+//orly5fj0aNHhVqH\nm5sbfH194ePjg3LlysHR0RFKpRIKhQIajQZqtRppaWmIj49HfHw84uLikJaWVqhtVa9eHaGhofjw\nww/h4eFRqHUwZqvS0tKwadMmhIWF4fLly4Vah6OjI3x8fHRt1sXFBUqlEkqlElqtFmq1GhkZGXj8\n+LGuvT5//rxQ2/Ly8sLw4cMxatQoVKpUqVDrYIwxa0JEiIyMRFhYGHbv3g2NRlPgdSgUCpQtW1bX\nF3t6esLOzg5KpRIAoFaroVKpkJSUhLi4OMTFxeHx48fQarUF3pZSqUSvXr0QGhqKtm3bmv0IY2FY\nbZKcQ61W4/Tp0wgPD0d4eDj+/fdfqUPSadCgAYKDgxEcHIxGjRpZZQUzZklEhGvXrmHv3r0IDw/H\n2bNn8zy8ZilvvvkmgoODERQUhLZt28Le3l7qkBhjzCyePHmCiIgIhIeH4+DBg3lOx7AUZ2dn3fS7\nbt26wdvbW+qQcmX1SbK+W7duYe/evdi3bx/Onz9v0Yp3d3dHs2bNEBQUhKCgIFSsWNFi22bMFiUk\nJOg66lOnTuHp06cW27aDgwMaNWqEbt26ITg4GHXr1uUfsoyxEicjIwORkZEIDw/H4cOHcfv2bYtu\nv1q1aujYsSOCg4PRrl07lCpVyqLbLwqbS5Jfp9FocOvWLcFcxOjoaKSmphZ53Z6enmjQoAEaN26M\nJk2aoHHjxqhSpQp/yTJWSESEhw8fCtprVFQUEhMTi7xuJycnwfy6xo0bo3bt2rCzszNB5IwxVnwk\nJycjKipK0BffuXOnUFMmXqdQKODv7y/ohxs2bGhVc4wLyqaTZGNSU1MF84zj4+ORmJgIlUoFlUoF\njUajm+tob2+PcuXKCeZC+vj4FOqkP8ZYwWVkZCAhIUHXVuPi4pCQkIDMzEyoVCqo1WooFAoolUrY\n2dnBy8tL0FZ9fX3h5ubGP2AZY6yQNBoNnjx5Isib4uLikJqaqpuHLJPJdLmTq6srfH19dX2xr68v\nvLy8it2Vgf4fnd40D7cnQvUAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"%%capture\n",
"pgm_markov_chain(\"images/pgm/markov-chain.png\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Markov Models: Joint Distribution\n",
"\n",
"If a sequence has the Markov property, then the joint distribution factorizes according to\n",
" $$\n",
" p(x_1, \\dots, x_N) = p(x_1) \\prod_{n=2}^N p(x_n | x_{n-1})\n",
" $$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Markov Models: Transition Matrix\n",
"\n",
"Suppose $X_t \\in \\{1,\\dots,K\\}$ is discrete. Then, a stationary chain with $N$ states can be described by a **transition matrix**, $A \\in \\R^{N \\times N}$ where\n",
" $$\n",
" a_{ij} = p(X_t=j \\mid X_{t-1}=i)\n",
" $$\n",
"\n",
"is the probability of transitioning from state $i$ to $j$.\n",
"\n",
"> Each row sums to one, $\\sum_{j=1}^K A_{ij} = 1$, so $A$ is a *row-stochastic matrix*."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Markov Models: Transition Diagram\n",
"\n",
"Transitions between states can be represented as a graph:\n",
"\n",
"\n",
"\n",
"> The nodes in this graph represent states that each $X_i$ can take.\n",
"\n",
"Here the transition matrix is $A = \\begin{bmatrix} (1-\\alpha) & \\alpha \\\\ \\beta & (1-\\beta) \\end{bmatrix}$\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Why do we care about Markov Models?\n",
"\n",
"**Answer**: they are the ideal model to reason about many randomized algorithms. Consider the following alg:\n",
"\n",
"```python\n",
"state = 0\n",
"for i in range(NUMITER):\n",
" cointoss = np.random.somedist() # sample a random variable\n",
" state = update(state,cointoss) # update state based on random sample\n",
"return state\n",
"```\n",
"\n",
"We can analyze the behavior of this algorithm by viewing the `state` variable as following a markov process.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Markov Models: State Vectors\n",
"\n",
"Consider a row vector $x_t \\in \\R^{1 \\times K}$ with entries $x_{tj} = p(X_t = j)$. Then,"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"$$\n",
" \\begin{align}\n",
" p(X_t = j)\n",
" &= \\sum_{i=1}^K p(X_t = j \\mid X_{t-1} = i) p(X_{t-1} = i) \\\\\n",
" &= \\sum_{i=1}^K A_{ij} x_{t-1,i} \\\\\n",
" \\end{align}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Therefore, we conclude $x_t = x_{t-1} A$.\n",
"\n",
"(Note $\\sum_{j=1}^K x_{tj} = 1$.)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Markov Models: Matrix Powers\n",
"\n",
"Since $x_t = x_{t-1} A$, this suggests that in general,\n",
" $$\n",
" x_{t} = x_{t-1} A = x_{t-2} A^2 = \\cdots x_{0} A^t\n",
" $$\n",
"If we know the initial state probabilities $x_0$, we can find the probabilities of landing in any state at time $t > 0$."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Example: Weather\n",
"\n",
"Suppose the weather is either $R=\\text{Rainy}$ or $S=\\text{Sunny}$,\n",
" $$\n",
" A = \\begin{bmatrix}\n",
" 0.9 & 0.1 \\\\\n",
" 0.5 & 0.5\n",
" \\end{bmatrix}\n",
" $$\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Example: Weather\n",
"\n",
"Suppose today is sunny, $x_0 = \\begin{bmatrix} 1 & 0 \\end{bmatrix}$. We can predict tomorrow's weather,\n",
" $$\n",
" x_1 = x_0 A = \\begin{bmatrix} 0.9 & 0.1 \\end{bmatrix}\n",
" $$\n",
"The weather over the next several days will be\n",
" $$\n",
" \\begin{align}\n",
" x_2 &= x_1 A = \\begin{bmatrix} 0.86 & 0.14 \\end{bmatrix} \\\\\n",
" x_3 &= x_2 A = \\begin{bmatrix} 0.844 & 0.156 \\end{bmatrix} \\\\\n",
" x_4 &= x_3 A = \\begin{bmatrix} 0.8376 & 0.1624 \\end{bmatrix}\n",
" \\end{align}\n",
" $$\n",
"\n",
"> **Question:** What happens to $x_0 A^n$ as $n \\rightarrow \\infty$?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Markov Chains: Stationary Distribution\n",
"\n",
"If we ever reach a stage $x$ where\n",
" $$\n",
" x = xA\n",
" $$\n",
"then we have reached the **stationary distribution** of the chain.\n",
"- To find $x=v^T$, solve the eigenvalue problem $A^T v = v$\n",
"- Under certain conditions, the limiting distribution $\\lim_{n \\rightarrow \\infty} x_0 A^n = x$\n",
"- Stationary distribution $x$ does not depend on the starting state $x_0$"
]
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}