teen_life.ipynb

gistfile1.txt

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src='tl.jpeg' width='600px' />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n_states = 5\n",
    "\n",
    "states = [\"Eat\", \"Hang\", \"Study\", \"Talk\", \"Sleep\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transition probability from state $i$ to state $j$  is $A_{i, j}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A = np.array([[0.,  0.7, 0.1, 0,   0.2],\n",
    "              [0.2, 0.4, 0,   0.2, 0.2],\n",
    "              [0,   1.0, 0,   0,   0],\n",
    "              [0,   0.3, 0,   0.7, 0],\n",
    "              [0.1, 0.1, 0,   0, 0.8]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.,  1.,  1.,  1.,  1.])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# probabilities must sum to one\n",
    "\n",
    "np.sum(A, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider an emsemble of teenagers. At each time $t$ each of them can be in either of the states *E*, *H*, *St*, *T*, *Sl*.\n",
    "\n",
    "A *population vector*, $s_{t}$, is a number of teenagers in each state at time $t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.068,  0.221,  0.155,  0.279,  0.277])"
      ]
     },
     "execution_count": 125,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.seed(1234)\n",
    "\n",
    "# Initial population vector\n",
    "s_0 = np.random.random(n_states)\n",
    "s_0 /= s_0.sum()\n",
    "\n",
    "s_0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time evolution of the population:\n",
    "\n",
    "$$\n",
    "s_{t+1} = A^T s_t\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.07186425,  0.40268538,  0.00679945,  0.23934926,  0.27930166])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s1 = np.dot(A.T, s)\n",
    "s1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "s_sleep = np.array([0, 0, 0, 0, 1.])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.1,  0.1,  0. ,  0. ,  0.8])"
      ]
     },
     "execution_count": 127,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(A.T, s_sleep)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import style as mstyle\n",
    "mstyle.use('ggplot')\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy import linalg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x5be28d0>"
      ]
     },
     "execution_count": 131,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEECAYAAAAmiP8hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1cVGXeP/DPzBmG4WEQB0VM0gyUxcmUQDKxtILbZdOV\n",
       "iviZZVu0r9aHfCjLQlMxIku41aw708K81+0Ju9Vq0/WFmmWSD6isCQbNqhkGIiALgjPDPPz+YDk4\n",
       "zgMD8nj4vF8vX84513Wduc4153w5c+ac85VZrVYriIhIsuRd3QEiIupYDPRERBLHQE9EJHEM9ERE\n",
       "EsdAT0QkcQz0REQS12Kgz8/Px4IFCzBv3jzs2LHDaT2dTodp06bh0KFDrW5LREQdx2Wgt1gsyMrK\n",
       "wuLFi7F69WocPHgQJSUlDut99NFHGD16dKvbXq+goKANqyFNHItmHItmHItmHAv3uAz0Op0OQUFB\n",
       "CAwMhEKhQExMDPLy8uzq7dq1C2PHjoWfn1+r216PH1wzjkUzjkUzjkUzjoV7XAb6qqoqBAQEiNMa\n",
       "jQZVVVV2dfLy8vBf//VfAACZTOZ2WyIi6ng3/GPs5s2bMX36dMhkMlitVvCJCkRE3YvCVaFGo0Fl\n",
       "ZaU4XVlZCY1GY1PnzJkzWLt2LQCgtrYW+fn5UCgUbrUFGr96Xfv1KykpqW1rIkEci2Yci2Yci2Yc\n",
       "C1vZ2dnia61WC61WC6CFQB8SEoKysjKUl5dDo9EgNzcX8+fPt6nzzjvviK/fffddREZGIioqCmaz\n",
       "ucW213emyW+//db6NZQgtVqN2traru5Gt8CxaMaxaMaxaHbTTTc5/cPnMtALgoDk5GSkp6fDYrHg\n",
       "vvvuQ3BwMHJycgAAcXFxrW5LRNSVLJYKmEzlDsuqjEZcunrZYVn/Om94VtY7LPP18IWPwsf+vQID\n",
       "YerXr+2dbSey7viYYh7RN+LRSjOORTOORbO2jIXRWIiaml0Oy85e9cI/zu11WHZv3Z3ok3PYYdmw\n",
       "vsPQ36u/3XxDfDyMI0a0qn9tddNNNzkt452xREQSx0BPRCRxLs/RE7UHtVp9w8sQBKFdliMF7TkW\n",
       "Uj4FVFGhQHm5/bGsWq1AXZ3gsE2DUtbR3eoSDPTUKaQcUHoqqf/hLC+XY9cuT7v5I0cqYDA4DvTq\n",
       "YR3dq67BUzdERBLHQE9EJHEM9EREEsdAT0QkcfwxljqdpcICU7mpw5avCFRA3s+9Y5jExEScOHEC\n",
       "gtD449zAgQPx7bfftvo9g4ODcfDgQQwZMqTVbYk6GgM9dTpTuQk1u2o6bPl+8X5Q9lO6XT89PR3T\n",
       "pk274ffthjeZEwHgqRsitwP02bNn8fDDDyM8PBwjR47E7NmzAQAPPfQQgMZnPw0fPhxfffUVACAn\n",
       "JwdxcXEYMWIEpk6ditOnT4vLuvPOO/HOO+/g3nvvhVarxfPPPw+DwQCgMZfDE088gREjRkCr1eKh\n",
       "hx7iHxG6IQz01OutXLkSI0eOREJCAn744Qen9TIyMjBx4kScPn0ax44dw1NPPQUA2LZtGwBgz549\n",
       "KC4uxpQpU3Dq1Cm88MILyMjIQEFBAR5//HE89dRTaGhoEJe3Y8cOfPzxxzh48CDOnDmDt956CwCw\n",
       "YcMG3HTTTfjxxx9x8uRJpKSkiAl9iNqCgZ56tSVLluDQoUM4fvw4HnvsMTz55JP45ZdfHNb18PDA\n",
       "r7/+itLSUiiVSowZM8bpcv/2t7/h8ccfx+jRoyGTyfDII49AqVTi+PHjABozsT355JMYOHAg/P39\n",
       "MW/ePHzxxRfi+5SXl+PXX3+FIAgu34fIHQz01KtFRETA29sbHh4eeOSRRzBmzBjs27fPYd1XXnkF\n",
       "ADB58mTcd999+Oyzz5wu98KFC9iwYQNGjBgh/istLUVZWZlY59qnDQ4aNEgsmzVrFm655RZMnz4d\n",
       "48aNw//8z/+0x6pSL8YfY4nc1L9/f6xatQoAcPToUUybNg1jx451eKXNTTfdhHnz5mHevHlOl3fh\n",
       "wgWb10FBQQAAHx8fLFu2DMuWLUNRURGSkpIwatQojB8/vp3XiHoLHtFTr1VTU4P9+/dDr9fDZDJh\n",
       "27ZtOHz4MCZOnOiw/ldffSXmSvDz8wMAyOWNu1D//v1tTvk89thj2LJlC06cOAGr1Yr6+nrs2bMH\n",
       "dXV1ABp/AP7f//1flJaW4vLly1i3bh3++Mc/Amj8Effs2bOwWq3w9fWFIAji5Z9ky1BqgLHQ6PCf\n",
       "+aIZ5kv2/6zG3vfDNo/oqdMpAhXwi/fr0OW7w2QyISMjAzqdDoIgIDQ0FJs2bcLQoUMd1j958iRS\n",
       "U1NRW1uL/v37Iy0tDTfffDMA4Pnnn8eCBQug1+uxatUqTJ48GRkZGXjllVdw9uxZqFQqREdHY9y4\n",
       "cQAaz9E/+OCDmD59Oi5evIhJkyaJqTbPnTuHpUuXorKyEn369MGf/vQn3HXXXe0wMtLj6lJdvZcM\n",
       "hp/t51ujrYBHB3esm2GGqW5MKpmEpLIe7Wns2LHIzMzs0tMxUvhc5GfkKN/uOC3gWa8++Mdeld38\n",
       "MY8VocHDcYYp9TAFDlza47CMGaaIiKjbavE7bn5+PjZv3iwm+E5ISLApP3r0KLKzsyGTySCTyTBj\n",
       "xgzcdtttAIA5c+bAy8sLcrkcgiBg5cqVHbMWRETklMtAb7FYkJWVhaVLl0Kj0SAlJQVRUVEIDg4W\n",
       "64wcOVK8zvf8+fPIyMjA22+/LZanpqbC19e3g7pP1DMdOnSoq7tAvYjLUzc6nQ5BQUEIDAyEQqFA\n",
       "TEwM8vLybOqoVM3nwPR6vXg1QpNu+BMAEVGv4vKIvqqqCgEBAeK0RqOBTqezq3fkyBF88sknuHz5\n",
       "snhTCdB4ZUFaWhrkcjliY2MRGxvbjl0nIiJ3tMvlldHR0YiOjsbp06fx9ttvi8/sSEtLQ9++fVFT\n",
       "U4O0tDQMGjQI4eHh7fGWRETkJpeBXqPRoLKyUpyurKyERqNxWj88PBwWiwW1tbVQq9Xo27cvgMab\n",
       "S6Kjo6HT6ewCfUFBAQoKCsTppKQkySctdpdSqZTEWPBmn+5JEIQev33pBT08lfYJwAFAISgcbnuC\n",
       "XIDFyTYpkzm/OU0Q5E7LFAoFlEr7R2PLVCp4duIYZ2dni6+1Wi20Wm1j/1w1CgkJQVlZGcrLy6HR\n",
       "aJCbmyve1NGkrKwMAwYMgEwmw5kzZwA0Xp9rMBhgsVjg5eUFvV6PkydPIjEx0e49ru1Mk55+bW97\n",
       "kcJ1zgB6fDCRKrPZ3OO3L7lZDoPR4LDMJKhgNtuHOLPFDLPZ7LCN1SpzWmY2W5yWmUwmGI1Gu/kG\n",
       "vR7GThpjtVqNpKQkh2UuA70gCEhOTkZ6erp4eWVwcDBycnIAND5/+/Dhw/juu+8gCAJUKpX4h6C6\n",
       "uhqZmZkAGq/eGT9+PEaNGtWe60VERG5o8Rx9REQEIiIibObFxcWJr6dOnYqpU6fatRswYAAyMjLa\n",
       "oYskNRUVCpSXd9y9eoGBFvTr516qwg8//BDZ2dkoKirC1KlTsWbNGpvyAwcOYMmSJfjtt98QERGB\n",
       "tWvXYtCgQa3qz2effYZPP/0U27dvb1U7ovbCZ91Qpysvl2PXLsfnVdtDfLwB/fq5VzcoKAgLFiwQ\n",
       "H252raqqKjzzzDPIzMxEXFwcVq1ahZkzZ4oZpIh6Cj4CgXq1+Ph4TJo0Sbxw4Fo7d+5EWFgYHnjg\n",
       "ASiVSixcuBCFhYX417/+5XBZn332GcaNG4ewsDDcdddd2L59O3Q6HVJSUnDs2DEMHz5c/D3KYDDg\n",
       "1VdfRXR0NEaPHo2XX35Z/EOTm5uLyMhIvP322xg5ciTGjh1r821g7969uPfeexEWFobIyEi89957\n",
       "HTAyJCU8oieC4xv7iouLMeKaB1J5eXlh6NChKCoqQkhIiE3d+vp6LF++HDt37sStt96KS5cu4fLl\n",
       "ywgNDcUbb7yBTz75xCZYv/766/j111+Rk5MDhUKBOXPmYM2aNUhJSQEAVFRU4PLlyzh+/DiOHTuG\n",
       "GTNmYNSoUbj11lvxwgsvYOPGjRgzZgxqampw/vz5DhoVkgoe0RMBDnOy1tfX210x5OvrKz5T/npy\n",
       "uRw//fQTrl69iv79+2P48OEA7P+IWK1WfPzxx1i+fDn69OkDHx8fPPvss/jyyy9t6r344ovw8PDA\n",
       "2LFjcf/994vlHh4eKCoqQm1tLfz8/MRnSxE5w0BPBMdH9N7e3naXH9bW1jp8dpO3tzfWr1+PLVu2\n",
       "IDIyEk888YTDu8iBxvtRrl69ivj4eDHN4IwZM1BVVSXW6dOnD7y8vMTp4OBgXLx4EQDw/vvvY9++\n",
       "fRg7diwSExNx7NixNq0z9R4M9ERwfEQfFhaGwsJCcbq+vh7nzp0Tj9SvN2HCBHzyySc4ceIEQkND\n",
       "sWjRIofL1mg0UKlU+Oabb1BYWIjCwkKcPn0aRUVFYp1///vfuHr1qjhdUlIiphocNWoUNm3ahJMn\n",
       "T2LSpEmYOXNm21ecegUGeurVzGYz9Ho9zObGm2gMBoN4U0x8fDyKioqwc+dO6PV6rF69Glqt1u78\n",
       "PNB4Tn337t2or6+Hh4cHvL29bdIMlpaWoqGhAUDjKZ7p06dj+fLl4p3npaWl+Pbbb22WmZmZiYaG\n",
       "Bhw6dAh79+7F5MmT0dDQgG3btqGmpgaCIIipBolc4Y+x1OkCAy2Ij3d8N2N7Ld9da9eutbl2ftu2\n",
       "bVi4cCGee+45aDQabNy4Ea+88grmzp2LO+64A+vXr3e4HIvFgvfffx8LFiyATCaDVqvFG2+8AQAY\n",
       "P348hg8fjtGjR0MQBJw8eRJLlizBmjVrMGXKFFRVVSEoKAh/+tOfMGHCBACNfxz69OmDO+64A97e\n",
       "3njzzTcREhIiBvqlS5fCbDYjNDTU5rHgRI4wlWA3JqVHIEhhPTpLbm4u5s2bZ/dI8PYmhc+FqQSb\n",
       "MZUgEVEvxkBP1A05+nGYqK0Y6Im6mXHjxuHo0aNd3Q2SEAZ6IiKJY6AnIpI4BnoiIoljoCcikjgG\n",
       "eiIiiWOgJyKSOD4CgTqdxVIBk8nx3YztQaEIhFzecoopo9GIlJQUfP/996iursaQIUOQkpKCe++9\n",
       "V6zDVIIkBS0G+vz8fGzevFlMDp6QkGBTfvToUWRnZ0Mmk0Emk2HGjBni87Fbaku9k8lUjpoax7eg\n",
       "twc/v3golS0HerPZjEGDBmHbtm0YNGgQ9uzZg5kzZ2Lv3r0IDg5mKkGSDJenbiwWC7KysrB48WKs\n",
       "Xr0aBw8eRElJiU2dkSNHIiMjA6tWrcKcOXOwYcMGt9uSa4ZSA4yFxlb/s1S4/1Cv3szLywvPP/+8\n",
       "eIQeGxuLm2++GT/++CMAphIk6XB5RK/T6RAUFITAwEAAQExMDPLy8hAcHCzWUamaHxqk1+vh5+fn\n",
       "dltyzVRuQs2umla384v3g7KfsgN6JG2XLl3C2bNnxefNM5UgSYXLI/qqqioEBASI0xqNxiYLTpMj\n",
       "R47gueeew+uvv46nnnqqVW2JuoOGhgY8++yzeOSRR8QgzlSCJBXt8mNsdHQ0oqOjcfr0abz99ttY\n",
       "u3at220LCgpQUFAgTiclJdntXL2VXtDDU+nZ6nYqlQo+ap8O6FHbdPfEGBaLBfPmzYOnpyfS09PF\n",
       "+W1JJfjee+/hhRdeQFRUFJYtW4bQ0FC7utemEmxitVphsTSfcmspleBbb72FlStXIjw8HCkpKYiM\n",
       "jGz1eguC0OP3NVf7iEJQONz2BLkAi5NtUiYTnG6vgiB3WqZQKKBU2n+LlqlU8OzEMc7OzhZfa7Va\n",
       "8VShy0Cv0WjEDDhA4waq0Wic1g8PD4fFYsGVK1fcbnttZ5r09Gdktxe5WQ6D0XGCjtoAH1TUO/74\n",
       "VL/IIOjt2wUGWtCvn6ld++iO7hxMrFYrFi5ciMrKSmzZssVmRw4LC8PWrVvFaXdSCU6YMAEGgwFv\n",
       "vvkmFi1ahG3btrlMJThgwACHy2pKJdgU7EtKShAeHg6gOZWg2WzGpk2bMHPmzDY9BM1sNvf4fc3V\n",
       "PmISVDCb7fcRs8UsZhG7ntUqc1pmNluclplMJhiNRrv5Br0exk4aY7VajaSkJIdlLgN9SEgIysrK\n",
       "UF5eDo1Gg9zcXMyfP9+mTllZGQYMGACZTIYzZ86Ib+jt7d1iW2q7inqFw6QKAHB7n2p4NVyymy/0\n",
       "qUPppQqHbfrXecOzst5hma+HL3wU9t8QLIGBMPVr+eqW7uzll1+GTqfDZ599Bk9P2yPD+Ph4vPba\n",
       "a9i5cyfuu+++FlMJHjt2DHfffTdUKpXTVIIeHh42qQTT09MREBCA0tJSFBcXixmmgMZUgi+//DKO\n",
       "HTuGvXv34sUXX0RDQwO++uorxMbGws/Pj6kEyS0uA70gCEhOTkZ6erp4iWRwcDBycnIAAHFxcTh8\n",
       "+DC+++47CIIAlUolBnNnbanjeXpdhMFgf/liaW3bs+d4OsmegzYEeoUiEH5+8S1XbCOFItCteiUl\n",
       "Jfjoo4/g6emJ0aNHi/NXrVqFhIQEphIkyWAqwW6sLWnSAOep0roqTZoUUtZ1JqYSdB9TCTZjKkEi\n",
       "ol6MgZ6oG2IqQWpPDPRE3QxTCVJ7Y6AnIpI4BnoiIoljoCcikjgGeiIiiWOgJyKSOAZ6IiKJYypB\n",
       "6nQV+gqU13dcKsFA70D0U7n3aIYPP/wQ2dnZKCoqwtSpU7FmzRqb8pZSCaanp+OTTz4BAEyfPh2L\n",
       "Fy9udX8TExPx8MMP49FHH211WyJ3MNBTpyuvL8eusx2XSjB+aLzbgT4oKAgLFizA/v37xQxPTVpK\n",
       "Jbhlyxbs3r0be/Y03jL/6KOP4uabb8aMGTPad4WIbhBP3VCvFh8fj0mTJqFv3752ZS2lEty6dStm\n",
       "zpyJoKAgBAUFYebMmTbPA7+WXq/H3Llzcdttt2HEiBF44IEHUFFRgTfeeANHjhzBK6+8guHDh2Pp\n",
       "0qUAGjO0TZs2DVqtFvfcc49NntoFCxbgpZdewqOPPoqwsDAkJibiwoULYvny5csxatQo/O53v0Ns\n",
       "bCyKiorac8ioB+IRPRHss0ABzlMJFhcXIyQkBD///LNNeXh4OIqLix0uf+vWraitrUVeXh48PT1R\n",
       "UFAAlUqFl19+GXl5eUhMTMS0adMAND73ftq0aVi0aBE+/vhjFBYW4tFHH8Xvfvc7DBs2DACwY8cO\n",
       "bNmyBaNHj8Zrr72GZ599Ftu3b8f+/ftx5MgRfP/991Cr1dDpdGJ6T+q9eERPBMfPlnGWSvDKlSsA\n",
       "gLq6OptyV2kGlUolLl++jLNnz0Imk+G2226zyVR17R+anJwcDB48GElJSZDL5bjtttsQHx+Pv//9\n",
       "72Kd2NhYREdHQ6lUis+sLy0thYeHB65cuYKff/4ZFosFoaGhYt5m6r14RE8Ex0f0LaUS9PHxEYN+\n",
       "U5mPj+MUjg8//DB+++03zJ49GzU1NXjooYfw0ksvQaFo3AWv/UNz4cIFnDhxwubbgslkQmJiolh3\n",
       "4MCBNv309/fHxYsXERMTg6eeegpLlixBSUkJ4uPjsWzZMofpD6n34BE9ERwf0YeFhaGwsFCcvj6V\n",
       "4PDhw23yHRcWFiIsLMzh8hUKBZ577jl88803+OKLL7Bnzx58/vnnDt970KBBGDt2LAoLC8V/xcXF\n",
       "eP311wE0/lG6NmdDXV0dqqurxbSEycnJ2LVrF/bv348zZ844TZZCvQcDPfVqZrMZer0eZnNjHlGD\n",
       "wSDmBY2Pj0dRURF27twJvV5vl0owMTERGzduRFlZGUpLS7Fx40anOTtzc3Nx+vRpmM1m+Pj4QKFQ\n",
       "2KQaPHfunFg3NjYWZ86cwf/93/+hoaEBDQ0NyM/Ph06nE+vs27cPR48ehdFoxKpVqxAZGYmBAwfi\n",
       "n//8J44fP46GhgZ4eXlBpVIx1SDx1A11vkDvQMQP7bhUgoHe7p+TXrt2rc2189u2bcPChQvx3HPP\n",
       "tZhKcMaMGTh//jzuv/9+AMBjjz2Gxx9/3OH7XLp0CS+//DJKS0vh4+ODP/7xj+KpmKeffhoLFizA\n",
       "li1b8PDDD+PVV1/Fxx9/jBUrVmDFihWwWCzQarVYvnw5gMZvAAkJCVi9ejWOHTuG22+/XUwnWFtb\n",
       "i9TUVJw/fx6enp6YOHEiZs2a1boBJMlhoKdO10/Vz+3r3DvawoULsXDhQqfld999N7799lun5UuW\n",
       "LMGSJUtafJ+pU6di6tSpDssiIyNx4MABm3khISH461//6nR5Go1GzEl7rfHjx4vX9RM1aTHQ5+fn\n",
       "Y/PmzWKC74SEBJvyAwcO4Msvv4TVaoWXlxf+/Oc/Y8iQIQCAOXPmwMvLC3K5HIIgYOXKlR2zFkS9\n",
       "SDdM80zdnMtAb7FYkJWVhaVLl0Kj0SAlJQVRUVEIDg4W6wwYMAArVqyAt7c38vPzsXHjRqSnp4vl\n",
       "qamp/MWfqB3JZDKmGqRWcRnodTodgoKCxOtwY2JikJeXZxPom65AAIDQ0FBUVlbaLINHH0Tt6/rn\n",
       "8RC1xGWgr6qqQkBAgDit0Whsfvm/3r59+xARESFOy2QypKWlQS6XIzY2FrGxse3QZSIiao12+zH2\n",
       "1KlT+Oabb5CWlibOS0tLQ9++fVFTU4O0tDQMGjQI4eHhNu0KCgpsrkVOSkqyuxuxt9ILengqPR2W\n",
       "KQSF08vmBLkAi4MymUxw3kaQOy1TKBRQKpX2y1Op4OnGZ8XL+7onQRB6/L7Wln3E2f4BdN0+0l6u\n",
       "fdaSVquFVqtt7J+rRhqNxuZUTGVlJTQajV29X375BRs2bMCSJUtszsc3PSjKz88P0dHR0Ol0doH+\n",
       "2s40uf5uxN5KbpbDYDQ4LDMJKpjNjj8+s8UsXgt+LatV5nA+AJjNFqdlJpMJRqPRbr5Br4fRjc+q\n",
       "pwcTqTKbzT1+X2vLPuJs/wC6bh9pD2q12ul9HC5vmAoJCUFZWRnKy8thMpmQm5uLqKgomzoVFRXI\n",
       "zMzE3LlzERQUJM43GAy4evUqgMYn9508eRKDBw++0XUhIqJWcnlELwgCkpOTkZ6eLl5eGRwcjJyc\n",
       "HABAXFwcPv/8c9TV1eGDDz4Q26xcuRLV1dXIzMwE0Hj1zvjx4zFq1KgOXh0iIrpei+foIyIibH5g\n",
       "BRoDfJOZM2di5syZdu0GDBiAjIyMdugiERHdCN4ZS51OUVEBeXnHpRK0BAbC1M+9O28TExNx4sQJ\n",
       "8Ue2gQMH2twJy1SCJAUM9NTp5OXl8NzVcakEDfHxgJuBHmgM1k1JP67FVIIkFXx6JfV6zm7qYypB\n",
       "kgoGeur1Vq5ciZEjRyIhIQE//PCDON9VKkEAbU4lWFBQgDfffFNMJRgdHY309HQUFxcjLS1NTCX4\n",
       "0EMP4ccff8S7776LxYsX4+effxaXt2PHDjz33HP48ccfMWLECDz77LMAYJNK8KeffsJ7773nMB8u\n",
       "9S4M9NSrLVmyBIcOHcLx48fx2GOP4cknn8T58+cBMJUgSQcDPfVqERER8Pb2hoeHBx555BGMGTMG\n",
       "e/fuBdD+qQQnTpyI2bNnIzIyEunp6TCZTGK5s1SCTf927NiBS5cuiXXdSSU4atQoLFq0yKaP1Dsx\n",
       "0BM5wVSCJBUM9NRr1dTUYP/+/dDr9TCZTNi2bRsOHz6MiRMnAmAqQZIOXl5Jnc4SGNh4CWQHLt8d\n",
       "JpMJGRkZ0Ol0EAQBoaGh2LRpE4YOHQoATCVIksFAT53O1K9fq65z7ygajQZff/21yzpMJUhSwFM3\n",
       "RD0Mk/lQazHQE/UwTCVIrcVTN0Q9DFMJUmvxiJ6ISOIY6ImIJI6BnohI4niOnjrFjeaNFQTBab7O\n",
       "Xucq0FDR0KomgkaA3IfHdb0VAz11uPZIQK1Wq3t8Iuv2Ij8jR/n21iVu8Yv3g3KEsoN6RN0d/8QT\n",
       "EUlci0f0+fn52Lx5s5gcPCEhwab8wIED+PLLL2G1WuHl5YU///nPGDJkiFttiYio47k8ordYLMjK\n",
       "ysLixYuxevVqHDx4ECUlJTZ1BgwYgBUrViAzMxMPP/wwNm7c6HZbIiLqeC6P6HU6HYKCgsTEBTEx\n",
       "McjLy0NwcLBYp+mRrQAQGhqKyspKt9sSUfuqDfBBRb39bq26qIIAx0+xDAy0oF8/k8MykgaXgb6q\n",
       "qgoBAQHitEajsXlU6vX27duHiIiINrUlohtXUa/AP/aq7OZ7nveE0N9xoI+PN3SHZ8xRB2q3q25O\n",
       "nTqFb775Bmlpaa1qV1BQYJO8ISkp6YYvxZMKvaCHp9LTYZlCUDh9zrggF2BxUCaTCc7bCHKnZQqF\n",
       "Akql/RUbMpUKnp30WSmVSm4X/9GW7UKhUMBD6eGwjUolg1rteHndXVvGwtn+AfTsfQSATXJ6rVYL\n",
       "rVbb2D9XjTQajXgqBgAqKyuh0Wjs6v3yyy/YsGEDlixZIqZZc7fttZ1pwsvoGsnNchiMBodlJkEF\n",
       "s9nxx2e2mB1ec261ypxei242W5yWmUwmGI1Gu/kGvR7GTvqseHlls7ZsFyHhF+DV/5LDNg2edcg9\n",
       "U+GwrH+dNzwr6x2W+Xr4wkdhnzrREhjY+CjqTtCWsXC2fwA9fx9xlvjGZaAPCQlBWVkZysvLodFo\n",
       "kJubi/nz59vUqaioQGZmJubOnYugoKBWtSWizuHpdREGwy6HZaW1Chy45PgZ9vfW3Yk+OYcdlg3r\n",
       "OwyeXv3t5hvi47tFvgFq5jLQC4KA5ORkpKeni5dIBgcHIycnBwAQFxeHzz//HHV1dfjggw/ENitX\n",
       "rnTaloiIOleL5+gjIiLEH1ibxMXFia9nzpyJmTNnut2WiIg6F++MJSKSOAZ6IiKJY6AnIpI4Bnoi\n",
       "IoljoCcikjgGeiIiiWOgJyKSOAZ6IiKJY6AnIpI4BnoiIoljoCcikjgGeiIiiWOgJyKSOAZ6IiKJ\n",
       "Y6AnIpI4BnoiIoljoCcikjgGeiIiiWOgJyKSuBZzxubn52Pz5s1igu+EhASb8gsXLuDdd9/FuXPn\n",
       "MG3aNEyZMkUsmzNnDry8vCCXy8Wk4URE1LlcBnqLxYKsrCwsXboUGo0GKSkpiIqKQnBwsFhHrVYj\n",
       "OTkZR48edbiM1NRU+Pr6tm+viYjIbS5P3eh0OgQFBSEwMBAKhQIxMTHIy8uzqePn54eQkBAIguBw\n",
       "GVartf16S0RErebyiL6qqgoBAQHitEajgU6nc3vhMpkMaWlpkMvliI2NRWxsbNt7SkREbdLiOfob\n",
       "kZaWhr4AQ+YXAAARyklEQVR9+6KmpgZpaWkYNGgQwsPDbeoUFBSgoKBAnE5KSoJare7IbvUYekEP\n",
       "T6WnwzKFoHD6LUqQC7A4KJPJBOdtBLnTMoVCAaVSab88lQqenfRZKZVKbhf/0Zbtwtk2AfTs7YJj\n",
       "YSs7O1t8rdVqodVqG/vnqpFGo0FlZaU4XVlZCY1G4/ab9u3bF0Dj6Z3o6GjodDq7QH9tZ5rU1ta6\n",
       "/R5SJjfLYTAaHJaZBBXMZscfn9lihtlstptvtcoczgcAs9nitMxkMsFoNNrNN+j1MHbSZ6VWq7ld\n",
       "/Edbtgtn2wTQs7cLjkUztVqNpKQkh2Uuz9GHhISgrKwM5eXlMJlMyM3NRVRUlMO615+LNxgMuHr1\n",
       "KgBAr9fj5MmTGDx4cFv6T0REN8DlEb0gCEhOTkZ6erp4eWVwcDBycnIAAHFxcaiurkZKSgrq6+sh\n",
       "l8uxc+dOrFmzBv/+97+RmZkJoPHqnfHjx2PUqFEdv0ZERGSjxXP0ERERiIiIsJkXFxcnvvb398f6\n",
       "9evt2qlUKmRkZLRDF4mI6EbwzlgiIoljoCcikjgGeiIiiWOgJyKSOAZ6IiKJY6AnIpI4BnoiIolj\n",
       "oCcikjgGeiIiiWOgJyKSOAZ6IiKJY6AnIpI4BnoiIoljoCcikjgGeiIiiWOgJyKSOAZ6IiKJY6An\n",
       "IpI4BnoiIolrMWdsfn4+Nm/eLCYHT0hIsCm/cOEC3n33XZw7dw7Tpk3DlClT3G5LREQdz+URvcVi\n",
       "QVZWFhYvXozVq1fj4MGDKCkpsamjVquRnJxsE+DdbUtERB3PZaDX6XQICgpCYGAgFAoFYmJikJeX\n",
       "Z1PHz88PISEhEASh1W2JiKjjuQz0VVVVCAgIEKc1Gg2qqqrcWvCNtCUiovbT4jn6jlZQUICCggJx\n",
       "OikpCWq1ugt71H3oBT08lZ4OyxSCwu5bVBNBLsDioEwmE5y3EeROyxQKBZRKpf3yVCp4dtJnpVQq\n",
       "uV38R1u2C2fbBNCztwuOha3s7GzxtVarhVarbeyfq0YajQaVlZXidGVlJTQajVtv6G7bazsj1j1c\n",
       "aVcPAKqMSly67Hig+wRXQqYqt5tv8TagxlLhsE3/Om94VtY7LPP18IWPwsd+eYGBMPXr57BNe5Ob\n",
       "5TAYDQ7LTIIKZrPjj89sMcNsNtvNt1plDucDgNlscVpmMplgNBrt5hv0ehhra511v12p1WrUdtJ7\n",
       "dXdt2S6cbRNAz94uOBbN1Go1kpKSHJa5DPQhISEoKytDeXk5NBoNcnNzMX/+fId1rVZrm9ter2ZX\n",
       "jcP5v3r1wT/2Og70Yx6rQIPHbrv56mEKHLi0x2Gbe+vuRJ+cww7LhvUdBk+v/nbzDfHxQCcFeiKi\n",
       "9uAy0AuCgOTkZKSnp4uXSAYHByMnJwcAEBcXh+rqaqSkpKC+vh5yuRw7d+7EmjVroFKpHLYlIqLO\n",
       "1eI5+oiICERERNjMi4uLE1/7+/tj/fr1brclIqLOxTtjiYgkjoGeiEjiGOiJiCSOgZ6ISOIY6ImI\n",
       "JI6BnohI4hjoiYgkjoGeiEjiGOiJiCSOgZ6ISOIY6ImIJI6BnohI4hjoiYgkjoGeiEjiGOiJiCSO\n",
       "gZ6ISOIY6ImIJI6BnohI4lpMJZifn4/NmzeLeV8TEhLs6mzatAn5+fnw9PTE7NmzMXToUADAnDlz\n",
       "4OXlBblcDkEQsHLlyvZfAyIicslloLdYLMjKysLSpUuh0WiQkpKCqKgomyTfx48fx8WLF7Fu3Tr8\n",
       "/PPP+OCDD5Ceni6Wp6amwtfXt+PWgIiIXHJ56kan0yEoKAiBgYFQKBSIiYlBXl6eTZ28vDxMmDAB\n",
       "ADBs2DDU1dWhurpaLLdarR3QbSIicpfLI/qqqioEBASI0xqNBjqdzmWdgIAAVFVVwd/fHzKZDGlp\n",
       "aZDL5YiNjUVsbGw7d5+IiFrS4jl6dzg7an/11Veh0WhQU1ODtLQ0DBo0COHh4e3xlkRE5CaXgV6j\n",
       "0aCyslKcrqyshEajcbtO0/9+fn6Ijo6GTqezC/QFBQUoKCgQp5OSkuCp9HTcWUEBQRAclglyARYH\n",
       "ZTKZ4LyNIHdaplAooFQq7ZenUsFTrXbYpr3pBT3H4j+USiXUnfRe3V1btgtn2wTQs7cLjoWt7Oxs\n",
       "8bVWq4VWq23sn6tGISEhKCsrQ3l5OTQaDXJzczF//nybOlFRUdi9ezdiYmJQXFwMHx8f+Pv7w2Aw\n",
       "wGKxwMvLC3q9HidPnkRiYqLde1zbmSYGo8Fhf0yCCmaz4y6bLWaYzWa7+VarzOF8ADCbLU7LTCYT\n",
       "jEaj3XyDXg9jba3DNu1NbpZzLP5DrVajtpPeq7try3bhbJsAevZ2wbFoplarkZSU5LDMZaAXBAHJ\n",
       "yclIT08XL68MDg5GTk4OACAuLg533HEHTpw4gblz50KlUmHWrFkAgOrqamRmZgJovHpn/PjxGDVq\n",
       "VHuuFxERuaHFc/QRERGIiIiwmRcXF2cz/fTTT9u1GzBgADIyMm6we0REdKN4ZywRkcQx0BMRSRwD\n",
       "PRGRxDHQExFJHAM9EZHEMdATEUkcAz0RkcQx0BMRSRwDPRGRxDHQExFJHAM9EZHEMdATEUkcAz0R\n",
       "kcQx0BMRSRwDPRGRxDHQExFJHAM9EZHEMdATEUlci6kEiboDQ6kBxvP2yZcBoMqoxKXLgsOyPsGV\n",
       "kKnK7eZbvA2osVQ4bNO/zhuelfUOy3w9fOGj8LFfXmAgTP36Oes+UZdqMdDn5+dj8+bNYnLwhIQE\n",
       "uzqbNm1Cfn4+PD09MXv2bAwdOtTttkTuMJWbULOrxmHZr1598I+9jgP9mMcq0OCx226+epgCBy7t\n",
       "cdjm3ro70SfnsMOyYX2HwdOrv918Q3w8wEBP3ZTLUzcWiwVZWVlYvHgxVq9ejYMHD6KkpMSmzvHj\n",
       "x3Hx4kWsW7cOzzzzDD744AO32xIRUcdzGeh1Oh2CgoIQGBgIhUKBmJgY5OXl2dTJy8vDhAkTAADD\n",
       "hg1DXV0dqqur3WpLREQdz2Wgr6qqQkBAgDit0WhQVVXlsk5AQACqqqrcaktERB2vXa66sVqt7bEY\n",
       "IiLqADKriyhdXFyMrVu3YsmSJQCA7du3QyaT2fyounHjRmi1WsTExAAAFixYgNTUVJSXl7fYFgAK\n",
       "CgpQUFAgTiclJbXf2hER9SLZ2dnia61WC61WC6CFI/qQkBCUlZWhvLwcJpMJubm5iIqKsqkTFRWF\n",
       "7777DkDjHwYfHx/4+/u71bapM0lJSeK/azva23EsmnEsmnEsmnEsbF0bS5uCPNDC5ZWCICA5ORnp\n",
       "6eniJZLBwcHIyckBAMTFxeGOO+7AiRMnMHfuXKhUKsyaNctlWyIi6lwtXkcfERGBiIgIm3lxcXE2\n",
       "008//bTbbYmIqHN1u0cgXPt1o7fjWDTjWDTjWDTjWLjH5Y+xRETU83W7I3oiImpfDPRERBLXJU+v\n",
       "/H//7/9hyJAh4nRMTAymTp3qtP7+/fsxatQo9O3btzO612FmzJiBLVu2iNP79+/HmTNnkJyc3IW9\n",
       "6nzbtm3DwYMHIZfLIZPJ8Mwzz6C4uBixsbFQKpWtWtb1Y9qT1dbWIi0tDQBQXV0NuVwOPz8/yGQy\n",
       "vP766xCE5ge3ZWdnw8vLC1OmTEFqaiqeeOIJ3HrrrV3V9XbjaNv429/+Jpn16ypdEug9PT2xatUq\n",
       "t+t/++23GDx4cI8P9DKZrKu70OWKi4tx/PhxvPnmm1AoFLhy5QqMRiO+/vpr3H333a0O9FIaU7Va\n",
       "Le4XW7duhZeXFyZPnuywrkwmE9ddKmPgaNtoaGiQzPp1pW71PPrPP/8cx44dg9FoRFhYGJ555hkc\n",
       "OnQI//rXv7Bu3Tp4enoiLS2t1cGgJ8jLy8P27dthMpng6+uLefPmoU+fPsjOzkZFRQUuXbqEiooK\n",
       "/OEPf0B8fDyAxvH6/vvv4efnh4CAANx6662YMmVKF6+Ja9XV1VCr1VAoGjc9X19f7Ny5E5cvX8aK\n",
       "FSvg5+eHZcuW2RypHzp0CMePH8fs2bNRXl6Ot956CwaDweYGvHfeeQd33nknxowZAwBYt24dxo0b\n",
       "5/AmvZ7CarVi79692LNnD0wmE4KCgjB37lyn27/FYsH69esREBCAadOmdXJvb5yjbeN6//znP7F1\n",
       "61Y0NDRgwIABmD17NlQqFc6cOYO//vWv0Ov1UKvVmDNnDvz9/ZGamopbbrkFhYWFMJvNmDVrFkJD\n",
       "Qzt71bpclwR6o9GIRYsWidMPPvgg7rrrLvz+979HYmIigMYd99ixYxg7dix2796NGTNm9Pivbtev\n",
       "95UrV8RAFB4eLr7eu3cvvvjiCzzxxBMAgNLSUixfvhxXr17FggULMGnSJJw5cwZHjhxBZmYmTCYT\n",
       "XnrpJYSEhHT+SrXS7bffjs8//xzz58/HyJEjMW7cOPzhD3/A119/jdTUVHHndnYU9+GHH2LSpEm4\n",
       "5557sHt383Pm77vvPnz99dcYM2YM6uvrUVxcjGeffbZT1qkj3Xnnnbj//vsBAJ9++in27duH3//+\n",
       "93b1zGYz1q1bhyFDhuDBBx/s7G62C0fbxogRI8TympoabNu2DcuWLYNSqcSOHTvw97//HQ8++CA2\n",
       "bdqEl156CWq1Grm5ufjkk08wa9YsyGQyGI1GrFq1CqdPn8b69evx3//93124ll2jSwK9Uql0eOrm\n",
       "1KlT+Oqrr2AwGHDlyhXcfPPNiIyM7IIedozr17vpHD0AVFZWYs2aNaiurobJZMKAAQMANAa8O+64\n",
       "AwqFAmq1Gn5+fqiurkZRURHGjBkDhUIBhUKByMjIHvFwOZVKhTfeeAM//fQTTp06hbVr12L69OkA\n",
       "3Hs4XnFxMV588UUAwN13342PPvoIADBixAhkZWWhpqYGhw4dwtixYyGX9/xrDc6fP49PP/0U9fX1\n",
       "0Ov1GD16tF0dq9WKjRs3Yty4cT02yAOutw0A+Pnnn1FSUoJXXnkFAGAymTB8+HD89ttv+PXXX8Xf\n",
       "NywWi81p3qbncIWHh+Pq1auor6+Ht7d3J65Z1+s2p26MRiOysrLw5ptvQqPRiF/PeotNmzZhypQp\n",
       "iIyMRGFhIbZu3SqWNX2VBQC5XA6z2QyZTGYTGHtCkG8il8sxYsQIjBgxAoMHD8b+/fsB2B7FX/va\n",
       "aHScQvB699xzD7777jv88MMPmD17drv2uau8++67WLRokThOhYWFdnVkMhnCwsJw6tQpTJ48GR4e\n",
       "Hl3Q0/Zx/bbx7bff2pTffvvtmD9/vs288+fP4+abb8Zrr73m1nv0xnP+3eaQpymo+/r6Qq/X49Ch\n",
       "Q2KZSqXC1atXu6prneLq1aviUUhT4AMcB/CmHfvYsWNoaGiAXq/H8ePHe8QG/Ntvv6G0tFScPnfu\n",
       "HAIDA+Hl5YX6+uY8rX369MGFCxdgsVhw5MgRcX5YWBgOHjwIAPj+++9tlj1x4kTs3LkTADBo0KCO\n",
       "XI1Oo9fr4e/vD5PJhAMHDojzrVarzbZx3333ISIiAmvWrIHFYumKrt4wR9tG//7NaRuHDRuGoqIi\n",
       "lJWVAWgcm9LSUtx0002oqalBcXExgMYj/Wuz2eXm5gIAfvrpJ/j4+MDLy6szVqdb6Rbn6EePHo3p\n",
       "06fj/vvvx8KFC+Hv72/zg8nEiRPx/vvv9/gfY10F4kceeQRr1qyBj48PbrvtNly6dEls46hdSEgI\n",
       "oqKi8MILL8Df3x+DBw/uEV9H9Xo9PvzwQ9TV1UEQBAQFBeEvf/kLvv/+e7z++uvQaDRYtmwZpk+f\n",
       "jjfeeAN+fn4ICQmBXq8HADz55JNYt24dvvjiC0RFRdmMTZ8+fRAcHIzo6OiuWr12l5SUhMWLF8PP\n",
       "zw+hoaHiODjaLiZPnoz6+nq8/fbbmDdvXo/4w38tR9vGM888g9WrVwMA/Pz8MHv2bLz11lswmUwA\n",
       "gGnTpmHgwIF4/vnn8eGHH6K+vh4WiwUPPPCA+BBFDw8PvPTSS+KPsb0RH4HQg+n1eqhUKhgMBqSm\n",
       "puIvf/kLbrnllq7uVpcxGAx44YUXsGrVql551Eb2VqxYIYkLOW5UtzlHT623ceNGlJSUoKGhARMm\n",
       "TOjVQf7kyZPYsGEDHnjgAQZ5ouvwiJ6ISOK6zY+xRETUMRjoiYgkjoGeiEjiGOiJiCSOgZ6ISOIY\n",
       "6ImIJO7/A4JMVJF147qsAAAAAElFTkSuQmCC\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x5e17690>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_steps = [5, 10, 20, 100, 500]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ccycle = ['m', 'b', 'y', 'g', 'r']\n",
    "\n",
    "for j, n in enumerate(n_steps):\n",
    "    \n",
    "    A_n = linalg.funm(A.T, lambda x: x**n)  # NB: this is *not*  (A.T)**n !\n",
    "    s = np.dot(A_n, s_0)\n",
    "\n",
    "    ax.bar(np.arange(n_states) + 0.1*j, s, \n",
    "           width=0.2, alpha=0.5, label=str(n)+' steps', color=ccycle[j])\n",
    "\n",
    "ax.set_xticks(np.arange(n_states))\n",
    "ax.set_xticklabels(actions)\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equilibrium distribuion:\n",
    "\n",
    "$$\n",
    "s_{e} = A^T s_e\n",
    "$$\n",
    "\n",
    "which corresponds to $A$ having at least one eigenvalue equal to $1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.9947164419982275e-18"
      ]
     },
     "execution_count": 132,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "linalg.det(A.T - np.eye(n_states))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.099,  0.297,  0.01 ,  0.198,  0.396])"
      ]
     },
     "execution_count": 206,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s_stable = linalg.solve(A.T - np.eye(n_states),\n",
    "                        np.ones(n_states))\n",
    "s_stable /= s_stable.sum()\n",
    "s_stable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.139+0.107j, -0.139-0.107j,  0.416+0.j   ,  1.000+0.j   ,\n",
       "        0.762+0.j   ])"
      ]
     },
     "execution_count": 133,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "linalg.eigvals(A.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Is the model ergodic? I.e. is there only one limiting distribution?|"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.099  0.099  0.099  0.099  0.099]\n",
      " [ 0.297  0.297  0.297  0.297  0.297]\n",
      " [ 0.01   0.01   0.01   0.01   0.01 ]\n",
      " [ 0.198  0.198  0.198  0.198  0.198]\n",
      " [ 0.396  0.396  0.396  0.396  0.396]]\n"
     ]
    }
   ],
   "source": [
    "B = A.T.copy()\n",
    "\n",
    "for _ in range(32):\n",
    "    B = np.dot(B, B)\n",
    "\n",
    "print B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "#### Simulate the random walk directly instead.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import scipy.stats as stats\n",
    "\n",
    "np.random.seed(1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 237,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# A[i,:] is a vector of transition probabilities from state `i` \n",
    "# prob[i] is a corresponding discrete probability distribution\n",
    "\n",
    "prob = [None] * n_states\n",
    "for i in range(n_states):\n",
    "    prob[i] = stats.rv_discrete(\n",
    "                    values=(np.arange(n_states), A[i, :])\n",
    "              )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 264,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.LineCollection at 0xb09cd10>"
      ]
     },
     "execution_count": 264,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEECAYAAADOJIhPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAHT9JREFUeJzt3X9MXfX9x/HX5ZIW2l4ClwSYkGqDNSNXU9FbbELVTaFZ\n",
       "9l0tRiVdp4lBo/1hbbPo1tJ2obtj3QRrtZVqUxr8w5ihicYluoVFk3Yl/kB7U3vXiHfVuuooKbd4\n",
       "EaET7vn+0XjDBXrvuT/ovfQ8H/94L/fzefO+t+f64t5zzufYDMMwBACwpKx0NwAASB9CAAAsjBAA\n",
       "AAsjBADAwggBALAwQgAALCw71gCv16uOjg6FQiHdcccdqqurm3ac3+/X9u3btXnzZi1btkyStGHD\n",
       "BuXm5iorK0t2u127du1KbffADGhra9OxY8eUl5enp59+etoxhw4dktfr1dy5c7V+/XotWrToMncJ\n",
       "pEbUTwKhUEjt7e1qbGzU7t27dfToUZ05c2bacS+//LJuvPHGKY81NTXpqaeeMh0APp/PZOvUotbM\n",
       "1PrpT3+qxsbGSz7+8ccf6+zZs3ruuef0yCOP6ODBgyn73WZRi1qpqhU1BPx+v0pKSlRUVKTs7GxV\n",
       "V1erp6dnyri3335by5YtU15e3pTH4j0XLRNeFGpZu1ZFRYXmz59/ycd7enp0++23S5IWL16s4eFh\n",
       "DQ4OpuR3m0UtaqWqVtQQCAQCKiwsDN93Op0KBAJTxvT09GjFihWSJJvNFn7MZrPJ4/Foy5Yt+sc/\n",
       "/pFQg0Cmmfy+KCwsnPK+AGaLmPsEYuno6NCaNWtks9lkGEbEX/4ej0cFBQUKBoPyeDwqLS1VRUVF\n",
       "sr8SSDtWW8GVwhZt7aDe3l69+uqr2rZtmyTp9ddfl81mi9g5/Nhjj4XfEENDQ5o7d64effRRud3u\n",
       "iFqvvvqqcnJytHLlyoif+3y+iI8x9fX1yT8rIIrOzs7wbZfLJZfLNWVMf3+//vznP0+7Y/jAgQNy\n",
       "uVyqrq6WJG3evFlNTU3Kz8+PGMe2jcvNzLY9WdRPAuXl5err61N/f7+cTqe6u7u1adOmiDH79u0L\n",
       "325ra9PNN98st9utCxcuKBQKKTc3V6Ojozp+/LjuvffeKb9juka//vrrmI2b4XA4NDQ0RC1qhV11\n",
       "1VVJ/8/Y7Xbr73//u6qrq9Xb26v58+dPCQCJbZtaM1urtLQ04r5hGAlt21FDwG63q6GhQc3NzeFD\n",
       "RMvKytTV1SVJqq2tveTcwcFBtba2Srp49NDy5cu1ZMmSuBsELrc9e/bo5MmTCgaDWrdune677z6N\n",
       "j49LurjN33TTTTp27Jg2btyonJwcrVu3Ls0dA4mL+nVQuvDXErVmqtZVV12VkjqJYtumVqpqTfdJ\n",
       "IBGcMQwAFkYIAICFEQIAYGGEAABYGCEAABZGCACAhRECAGBhhAAAWBghAAAWRggAgIURAgBgYYQA\n",
       "AFgYIQAAFkYIAICFEQIAYGGEAABYWMwQ8Hq92rx5sx5//HG98cYblxzn9/u1evVqvffee3HPBQCk\n",
       "R9QQCIVCam9vV2Njo3bv3q2jR4/qzJkz0457+eWXdeONN8Y9FwCQPlFDwO/3q6SkREVFRcrOzlZ1\n",
       "dbV6enqmjHv77be1bNky5eXlxT0XAJA+UUMgEAiosLAwfN/pdCoQCEwZ09PToxUrVkiSbDab6bkA\n",
       "gPTKTrZAR0eH1qxZI5vNJsMw4r7Ysc/nk8/nC9+vr6+Xw+FIti1J0pw5cxKuNfFTzXSCwWBCdaXk\n",
       "+qJW8jo7O8O3XS6XXC5XymoDs03UEHA6nRoYGAjfHxgYkNPpjBhz6tQp7dmzR5I0NDQkr9er7Oxs\n",
       "U3Ol6d+EQ0ND8T+TaTgcjpTVmiyZuqnsi1rx16qvr09JLeBKEDUEysvL1dfXp/7+fjmdTnV3d2vT\n",
       "pk0RY/bt2xe+3dbWpptvvllut1vj4+Mx5wIA0itqCNjtdjU0NKi5uVmhUEh33HGHysrK1NXVJUmq\n",
       "ra2Ney4AIHPE3CdQWVmpysrKiJ9d6n/+69evjzkXAJA5OGMYACyMEAAACyMEAMDCCAEAsDBCAAAs\n",
       "jBAAAAsjBADAwggBALAwQgAALIwQAAALIwQAwMKSvp4AcKXxer3q6OgIL3xYV1cX8XgwGNTevXs1\n",
       "ODioUCiklStX6ic/+Ul6mgWSRAgAE/xwbewdO3bI6XRq69atcrvdESvg/u1vf9OiRYu0Zs0aBYNB\n",
       "bd68WbfeeqvsdnsaOwcSw9dBwARmro1dUFCg7777TpI0MjIih8NBAGDWIgSACcxcG/vOO+/UmTNn\n",
       "9Oijj+rJJ5/Ugw8+eJm7BFKHr4OAOL3++uu65ppr1NTUpL6+Pv3hD39QS0uLcnNzI8Zl6vWzqRWf\n",
       "2XK9cSmx62fHDIFYO8k+/PBDdXZ2ymazyWaz6YEHHtD1118vSdqwYYNyc3OVlZUlu92uXbt2xfuc\n",
       "gMvKzLWxe3t7dffdd0tS+Kujr7/+WuXl5RHjZsv1s6mVnEy53rikhK6fHTUEzOwku+GGG7R06VJJ\n",
       "0pdffqmWlhbt3bs3/HhTU5MWLFgQd2NAOpi5rvZVV12lTz75RD/+8Y81ODior7/+WsXFxWnqGEhO\n",
       "1BCYuJNMUngn2cQQyMnJCd8eHR2d8tHJMIxU9gvMKDPX1b777rvV1tamJ598UqFQSPfffz9/6GDW\n",
       "ihoC0+0k8/v9U8Z98MEHeuWVV3T+/Hlt3749/HObzSaPx6OsrCzV1NSopqYmha0DMyPWdbXz8vK0\n",
       "ZcuWy90WMCNSsmO4qqpKVVVVOnnypPbu3atnn31WkuTxeFRQUKBgMCiPx6PS0lJVVFREzJ0tO88m\n",
       "S6ZupuzwsmItKbGdZ8CVKmoImNlJNlFFRYVCoZCGhobkcDhUUFAg6eJfTlVVVfL7/VNCYLbsPJss\n",
       "U3YGUSv+WonsPAOuVFHPE5i4k2xsbEzd3d1yu90RY/r6+sLf+586dUrSxTfahQsXNDIyIunivoLj\n",
       "x49r4cKFM/EcAAAJivpJwMxOsvfff1+HDx+W3W5XTk5O+EiKwcFBtba2Srp4lNHy5cu1ZMmSGX46\n",
       "AIB4xNwnEGsn2apVq7Rq1aop84qLi9XS0pKCFgEAM4VlIwDAwggBALAwQgAALIwQAAALIwQAwMII\n",
       "AQCwMEIAACyMEAAACyMEAMDCCAEAsDBCAAAsjBAAAAsjBADAwggBALAwQgAALIwQAAALi3lRGa/X\n",
       "q46OjvCVxerq6iIe//DDD9XZ2SmbzSabzaYHHnhA119/vam5AID0ihoCoVBI7e3t2rFjh5xOp7Zu\n",
       "3Sq3262ysrLwmBtuuEFLly6VJH355ZdqaWnR3r17Tc0FAKRX1K+D/H6/SkpKVFRUpOzsbFVXV6un\n",
       "pydiTE5OTvj26Oio8vLyTM8FAKRX1E8CgUBAhYWF4ftOp1N+v3/KuA8++ECvvPKKzp8/r+3bt8c1\n",
       "FwCQPjH3CZhRVVWlqqoqnTx5Unv37tWePXtMz/X5fPL5fOH79fX1cjgcqWhLc+bMSVmtyZKpm8q+\n",
       "qBW/zs7O8G2XyyWXy5Wy2sBsEzUEnE6nBgYGwvcHBgbkdDovOb6iokKhUEjffvut6bnTvQmHhoZM\n",
       "P4FoHA5HympNlkzdVPZFrfhr1dfXp6QWcCWIuk+gvLxcfX196u/v19jYmLq7u+V2uyPG9PX1yTAM\n",
       "SdKpU6ckXXyjmZkLAEivqJ8E7Ha7Ghoa1NzcHD7Ms6ysTF1dXZKk2tpavf/++zp8+LDsdrtycnK0\n",
       "adOmqHMBAJkj5j6ByspKVVZWRvystrY2fHvVqlVatWqV6bkAgMzBGcMAYGGEAABYGCEAABZGCACA\n",
       "hRECAGBhKTljGLiSmFn91ufz6aWXXtL4+LgcDoeampouf6NAChACwARmVr8dHh5We3u7tm3bpsLC\n",
       "QgWDwTR2DCSHr4OACcysfvvPf/5Tt9xyS3iBxB9WzgVmIz4JABOYWf32v//9r8bHx7Vz506NjIzo\n",
       "5z//uW677bbL3SqQEoQAEKfx8XF9/vnn+t3vfqcLFy5o+/btWrx4sX70ox+luzUgboQAMIGZ1W8L\n",
       "CwvlcDg0Z84czZkzRxUVFTp9+vSUEMjUZdJjfX2VzD6OK7GvWDJlWXkpsWXSCQFggomr3zqdTnV3\n",
       "d4cXRfzB0qVLdejQIYVCIX3//ff67LPP9Itf/GJKLZZJTy36ii2RZdIJAWACMyvnlpaWasmSJXri\n",
       "iSdks9l05513skIuZi1CAJgk1sq5knTXXXfprrvuupxtATOCQ0QBwMJifhKIdfbkkSNH9Oabb8ow\n",
       "DOXm5urhhx/W1VdfLUnasGGDcnNzlZWVJbvdrl27ds3MswAAJCRqCJg5e7K4uFg7d+7UvHnz5PV6\n",
       "deDAATU3N4cfb2pq0oIFC2buGQAAEhb16yAzZ09ed911mjdvniTp2muvjTi8TlL4+sMAgMwT9ZOA\n",
       "mbMnJ3rnnXcidqjZbDZ5PB5lZWWppqZGNTU1KWgZAJAqKTs66MSJE3r33Xfl8XjCP/N4PCooKFAw\n",
       "GJTH41FpaakqKipS9SsBAEmKGgJmzp6UpNOnT+vFF1/Utm3bIr7/LygokHTxTMCqqir5/f4pIZCp\n",
       "Z1XGkilnCVIrfomcVQlcqaKGgJmzJ8+dO6fW1lZt3LhRJSUl4Z9fuHBBoVBIubm5Gh0d1fHjx3Xv\n",
       "vfdO+R2cVZkcasVfK5GzKoErVdQQMHP25Guvvabh4WEdPHgwPGfXrl0aHBxUa2urpItHGS1fvlxL\n",
       "liyZ4acDAIhHzH0Csc6eXLt2rdauXTtlXnFxsVpaWlLQIgBgpnDGMABYGCEAABZGCACAhRECAGBh\n",
       "hAAAWBghAAAWRggAgIURAgBgYYQAAFgYIQAAFkYIAICFEQIAYGGEAABYGCEAABZGCACAhRECAGBh\n",
       "MS8q4/V61dHREb6yWF1dXcTjR44c0ZtvvinDMJSbm6uHH35YV199tam5AID0ivpJIBQKqb29XY2N\n",
       "jdq9e7eOHj2qM2fORIwpLi7Wzp071draqnvuuUcHDhwwPRcAkF5RQ8Dv96ukpERFRUXKzs5WdXW1\n",
       "enp6IsZcd911mjdvniTp2muv1cDAgOm5AID0ihoCgUBAhYWF4ftOp1OBQOCS4995553w9YjjnQsA\n",
       "uPxStmP4xIkTevfdd/WrX/0qVSUBADMs6o5hp9MZ/npHkgYGBuR0OqeMO336tF588UVt27ZNCxYs\n",
       "iGuuz+eTz+cL36+vr1dpaWn4fjAYjOPpRJozZ44cDkfC86NJpm4q+8qUWnl5eVEfT9e/4+S+DMNQ\n",
       "Z2dn+L7L5ZLL5Uq4N2C2ixoC5eXl6uvrU39/v5xOp7q7u7Vp06aIMefOnVNra6s2btyokpKSuOZK\n",
       "sd+EQ0ND8T6nMIfDkdT8aDKlr0ytNVmmvF7SxT80AFwUNQTsdrsaGhrU3NwcPsyzrKxMXV1dkqTa\n",
       "2lq99tprGh4e1sGDB8Nzdu3adcm5AIDMEfM8gcrKyvDO3h/U1taGb69du1Zr1641PRcAkDk4YxgA\n",
       "LIwQACbxer3avHmzHn/8cb3xxhuXHOf3+7V69Wq9//77l7E7ILUIAWACs2e6h0Ihvfzyy7rxxhtl\n",
       "GEYaOgVSgxAAJjB7pvvbb7+tZcuWxTw0Fsh0hAAwgZkz3QOBgHp6erRixQpJks1mu6w9AqlECABx\n",
       "6ujo0Jo1a2Sz2WQYBl8HYVaLeYgoYCVmznQ/deqU9uzZI+niSXBer1fZ2dlyu90R46Y7Gz4Tzu6O\n",
       "JVPOhp+MvmJL5Gx4QgCYwMyZ7vv27Qvfbmtr08033zwlAKTp34Sc3Z04+ootkbPhCQFgAjNnyQNX\n",
       "EkIAmCTWWfITrV+//nK0BMwYdgwDgIURAgBgYYQAAFgYIQAAFkYIAICFEQIAYGExDxH1er3q6OgI\n",
       "HzNdV1cX8fhXX32ltrY2ffHFF1q9erVWrlwZfmzDhg3Kzc1VVlZW+IpjAIDMETUEflhWd8eOHXI6\n",
       "ndq6davcbnfEZSIdDocaGhr04YcfTlujqakpfPF5AEBmifp1kJlldfPy8lReXi673T5tDRbXAoDM\n",
       "FfWTwHTL6vr9ftPFbTabPB6PsrKyVFNTo5qamsQ7BQCk3IwuG+HxeFRQUKBgMCiPx6PS0lJVVFRE\n",
       "jJlupcWJMmmFvokypa9MrTVZprxeUmIrLQJXqqghYGZZ3WgKCgokXfzKqKqqSn6/f0oIxHoTZtIK\n",
       "fRNlSl+ZWmuyTHm9pMRWWgSuVFH3CUxcVndsbEzd3d3TLpkrTf3u/8KFCxoZGZEkjY6O6vjx41q4\n",
       "cGGK2gYApELUTwJmltUdHBzU1q1b9d133ykrK0tvvfWWnnnmGX3zzTdqbW2VdPEoo+XLl2vJkiUz\n",
       "/4wAAKbF3CcQa1nd/Px87d+/f8q8nJwctbS0pKBFAMBM4YxhALAwQgAALIwQAAALIwQAwMIIAQCw\n",
       "MEIAACyMEAAACyMEAMDCCAEAsDBCAAAsjBAAAAsjBADAwggBALAwQgAALIwQAAALIwQAwMJiXlTG\n",
       "6/Wqo6MjfGWxurq6iMe/+uortbW16YsvvtDq1au1cuVK03MBAOkV9ZNAKBRSe3u7GhsbtXv3bh09\n",
       "elRnzpyJGONwONTQ0BDxP3+zcwEA6RU1BPx+v0pKSlRUVKTs7GxVV1erp6cnYkxeXp7Ky8tlt9vj\n",
       "ngsASK+oIRAIBFRYWBi+73Q6FQgETBVOZi4A4PKIuU9gpvl8Pvl8vvD9+vr6iMcdDkfCtefMmZPU\n",
       "/Ggypa9MrTVZprxektTZ2Rm+7XK55HK5UlYbmG2ihoDT6dTAwED4/sDAgJxOp6nCZufGehMODQ2Z\n",
       "+n3TcTgcSc2PJlP6ytRak2XK6yVN/UNjslgHNBw5ckRvvvmmDMNQbm6uHn74YV199dUp6w+4nKJ+\n",
       "HVReXq6+vj719/drbGxM3d3dcrvd0441DCPhuUCmMHNAQ3FxsXbu3KnW1lbdc889OnDgQJq6BZIX\n",
       "9ZOA3W5XQ0ODmpubw38VlZWVqaurS5JUW1urwcFBbd26Vd99952ysrL01ltv6ZlnnlFOTs60c4FM\n",
       "NvGABknhAxombrvXXXdd+Pa1114b8YkXmG1i7hOorKxUZWVlxM9qa2vDt/Pz87V//37Tc4FMNt0B\n",
       "DX6//5Lj33nnHbZxzGpp3zEMzFYnTpzQu+++K4/HM+3j0x30UFpaGr4fDAYT/t1W2Yk/EX3FlshB\n",
       "D4QAMIHZAxpOnz6tF198Udu2bdOCBQumrcVBD6lFX7HFOuhhOqwdBExg5oCGc+fOqbW1VRs3blRJ\n",
       "SUmaOgVSg08CwARmDoZ47bXXNDw8rIMHD4bn7Nq1K51tAwkjBIBJYh0MsXbtWq1du/ZytwXMCL4O\n",
       "AgALIwQAwMIIAQCwMEIAACyMEAAACyMEAMDCCAEAsDBCAAAsjBAAAAsjBADAwggBALCwmGsHxbre\n",
       "qiQdOnRIXq9Xc+fO1fr167Vo0SJJ0oYNG5Sbm6usrCwW2QKADBQ1BH643uqOHTvkdDq1detWud3u\n",
       "iEvtffzxxzp79qyee+45ffbZZzp48KCam5vDjzc1NV1yvXUAQHpF/Tpo4vVWs7Ozw9dbnainp0e3\n",
       "3367JGnx4sUaHh7W4OBg+PHJF6AHAGSOqJ8EzFxvdfKYwsJCBQIB5efny2azyePxKCsrSzU1Naqp\n",
       "qUlx+wCAZKTkegKX+mv/97//vZxOp4LBoDwej0pLS1VRURExZrrrsE6USdfvnChT+srUWpNlyusl\n",
       "JXYdVuBKFTUEzFxvNdqYH/6bl5enqqoq+f3+KSHAdViTk6m1JsuU10tK7DqswJUq6j4BM9dbdbvd\n",
       "Onz4sCSpt7dX8+fPV35+vi5cuKCRkRFJ0ujoqI4fP66FCxfO0NMAACQi6icBM9dbvemmm3Ts2DFt\n",
       "3LhROTk5WrdunSRpcHBQra2tki4eZbR8+XItWbJkhp8OACAeMfcJxLreqiQ99NBDU+YVFxerpaUl\n",
       "yfYAADOJM4YBwMIIAQCwMEIAACyMEAAACyMEAMDCCAEAsDBCAAAsjBAAAAsjBADAwggBALAwQgAA\n",
       "LIwQAAALIwQAwMIIAQCwMEIAACws5vUEvF6vOjo6wheVqaurmzLm0KFD8nq9mjt3rtavX69FixaZ\n",
       "ngtkmmS2eWC2ifpJIBQKqb29XY2Njdq9e7eOHj2qM2fORIz5+OOPdfbsWT333HN65JFHdPDgQdNz\n",
       "gUyTzDYPzEZRQ8Dv96ukpERFRUXKzs5WdXW1enp6Isb09PTo9ttvlyQtXrxYw8PDGhwcNDUXyDTJ\n",
       "bPPAbBQ1BAKBgAoLC8P3nU6nAoFA1DGFhYUKBAKm5gKZJpltHpiNUrJj2DCMVJQBZg22eVwpou4Y\n",
       "djqdGhgYCN8fGBiQ0+k0NWZsbCzmXEny+Xzy+Xzh+/X19Sl9gzkcjoTmzfSbPNG+MrVWpr5e0/XV\n",
       "2dkZvu1yueRyucL3k9nmJ2Pbjg99xSfebTtaoUsaGxszHnvsMePs2bPG999/bzzxxBPGf/7zn4gx\n",
       "H330kfHHP/7RMAzD+PTTT43GxkbTc6fzl7/8JeYYs6hFrXhrJbPNJ/u740EtaqWqVtRPAna7XQ0N\n",
       "DWpubg4fLldWVqauri5JUm1trW666SYdO3ZMGzduVE5OjtatWxd1LpDJktnmgdko5nkClZWVqqys\n",
       "jPhZbW1txP2HHnrI9Fwg0yWzzQOzjb2pqakp3U1MVlRURC1qzYpa6fzd1KJWKmrZDIPDHADAqlg7\n",
       "CAAsjBAAAAuLuWN4pqRyka5YtXw+n5566ikVFxdLkm655Rbdc889U+q0tbXp2LFjysvL09NPPz3t\n",
       "7zLbU6xaZnuSpHPnzun555/XN998I5vNpjvvvFM///nPE+rNTC2zvf3vf/9TU1OTvv/+e42NjWnp\n",
       "0qVas2ZNQn2ZqRXPayZdXAdoy5Ytcjqd2rJlS0J9JYJtm2073lpp3bZTdpBqHMbHx+M6Fru3t/eS\n",
       "x2KbqXXixAnjT3/6U8y+/vWvfxmnTp0yfv3rX0/7uNmezNQy25NhGMb58+eNzz//3DAMwxgZGTEe\n",
       "f/zxhF8vM7Xi6W10dNQwjIvH1zc2NhonT55MqC8zteLpyzAM469//avx7LPPTjsnnr7iwbbNtp1I\n",
       "rXRu22n5OiiVi3SZXajOMLH/u6KiQvPnz7/k4/EsHBarltmeJCk/P1/XXHONJCknJ0elpaU6f/58\n",
       "Qr2ZqRVPb3PnzpUkjY2NKRQKacGCBQn1ZaZWPH0NDAzo2LFjuuOOO6adM1OLwLFtm+9JYttOpK9U\n",
       "b9tp+TpoukW6/H5/1DE/LNKVn58fdy2bzabe3l49+eSTcjqdeuCBBxI6cc1sT2Yk2lN/f7+++OIL\n",
       "LV68OOneLlUrnt5CoZB++9vf6uzZs1qxYsWUcfH0FatWPH299NJLuv/++zUyMjLt46n8t4xWl22b\n",
       "bdtMrXRu2xm9Y9hsMsayaNEi7d+/Xy0tLfrZz36mlpaWWdnT6Oiodu/erQcffFA5OTlJ9RatVjy9\n",
       "ZWVlqaWlRS+88IJOnjwZsVZOvH3FqmW2r48++kh5eXlatGhR1N+dqn/LRLBtR2LbTt+2nZYQSOUi\n",
       "XWbG5ebmhj+OVVZWamxsTN9+++2M9G1WvD2NjY3p6aef1q233qqqqqqkeotVK5HXa968eaqsrNS/\n",
       "//3vhPuKVctsX59++qk++ugjbdiwQc8++6x8Pp/27duXdF9msG2zbSdSK53bdlpCoLy8XH19ferv\n",
       "79fY2Ji6u7vldrsjxrjdbh0+fFiS1Nvbq/nz50/7ccZMrcHBwXAy/vBxerrv5GIx25MZ8fRkGIZe\n",
       "eOEFlZaW6v/+7/+S6s1MLbO9BYNBDQ8PS7p4BMQnn3wy5SgEs32ZqWW2rzVr1mj//v16/vnntXnz\n",
       "ZrlcLj322GMJ9RUvtm227URqpXPbTss+gVQu0mWm1nvvvaeuri5lZWVp7ty52rRp07S19uzZo5Mn\n",
       "TyoYDGrdunW67777ND4+HndPZmqZ7Um6mP5HjhzRwoUL9Zvf/EaS9Mtf/lLnzp2Luzcztcz2Njg4\n",
       "qOeff16hUEiGYei2227TDTfckNC/o5la8bxmE9lsNkm6LIvAsW2zbSdSK53bNstGAICFZfSOYQDA\n",
       "zCIEAMDCCAEAsDBCAAAsjBAAAAsjBADAwggBALAwQgAALOz/Aa1KH6aDqRfiAAAAAElFTkSuQmCC\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0xaab18d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2)\n",
    "\n",
    "x = np.arange(n_states)\n",
    "ax1.vlines(x, 0, prob[1].pmf(x), lw=8)\n",
    "ax2.vlines(x, 0, prob[1].cdf(x), lw=8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 269,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "num_steps = 10000\n",
    "histogram = np.zeros(n_states)\n",
    "\n",
    "state = 0\n",
    "for t in range(num_steps):\n",
    "    state = dists[state].rvs()\n",
    "    \n",
    "    histogram[state] += 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 270,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.099,  0.295,  0.012,  0.203,  0.392])"
      ]
     },
     "execution_count": 270,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "histogram / histogram.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.099,  0.297,  0.01 ,  0.198,  0.396])"
      ]
     },
     "execution_count": 213,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s_stable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "BUT: Drawing random variates one by one is *slow*."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}