{ "metadata": { "name": "", "signature": "sha256:8d9eb05a12220c324b7369687acad8b7fb9ce4a9960812c0537127537037ed80" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Exercises in Jarvis and Shier" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Reading:\n", "J. P. Jarvis and D. R. Shier, \"[Graph-Theoretic Analysis of Finite Markov Chain](http://www.ces.clemson.edu/~shierd/Shier/markov.pdf).\"" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from __future__ import division\n", "import numpy as np\n", "import quantecon as qe" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 6 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make sure that you have installed `quantecon` version 0.1.6 (or above):" ] }, { "cell_type": "code", "collapsed": false, "input": [ "qe.__version__" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 7, "text": [ "'0.1.6'" ] } ], "prompt_number": 7 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Exercise 3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider the Markov chain given by the following stochastic matrix\n", "(where the actual values of non-zero probabilities are not important):" ] }, { "cell_type": "code", "collapsed": false, "input": [ "P = np.zeros((6, 6))\n", "P[0, 0] = 1\n", "P[1, 4] = 1\n", "P[2, [2, 3, 4]] = 1/3\n", "P[3, [0, 5]] = 1/2\n", "P[4, [1, 4]] = 1/2\n", "P[5, [0, 3]] = 1/2" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 52 }, { "cell_type": "code", "collapsed": false, "input": [ "print (P)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 1. 0. ]\n", " [ 0. 0. 0.333 0.333 0.333 0. ]\n", " [ 0.5 0. 0. 0. 0. 0.5 ]\n", " [ 0. 0.5 0. 0. 0.5 0. ]\n", " [ 0.5 0. 0. 0.5 0. 0. ]]\n" ] } ], "prompt_number": 55 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a `MarkovChain` instance:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "mc0 = qe.MarkovChain(P)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 9 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We call the states $0, \\ldots, 5$, respectively, instead of $1, \\ldots, 6$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(a) Determine the communication classes." ] }, { "cell_type": "code", "collapsed": false, "input": [ "mc0.communication_classes" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 10, "text": [ "[array([0], dtype=int32),\n", " array([1, 4], dtype=int32),\n", " array([3, 5], dtype=int32),\n", " array([2], dtype=int32)]" ] } ], "prompt_number": 10 }, { "cell_type": "markdown", "metadata": {}, "source": [ "(b) Classify the states of this Markov chain." ] }, { "cell_type": "code", "collapsed": false, "input": [ "mc0.recurrent_classes" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 16, "text": [ "[array([0], dtype=int32), array([1, 4], dtype=int32)]" ] } ], "prompt_number": 16 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obtain a list of the recurrent states." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Write your own code\n", "RS=np.array((0, 1, 4))\n", "print(RS)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[0 1 4]\n" ] } ], "prompt_number": 44 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obtain a list of the transient states." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Write your own code\n", "TS=np.array((2, 3, 5))\n", "print(TS)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[2 3 5]\n" ] } ], "prompt_number": 45 }, { "cell_type": "markdown", "metadata": {}, "source": [ "(c) Does the chain have a unique stationary distribution?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "mc0.num_recurrent_classes" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 17, "text": [ "2" ] } ], "prompt_number": 17 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obtain the stationary distributions:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "vecs = mc0.stationary_distributions\n", "print (vecs)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0.33333333 0. 0. 0.66666667 0. ]]\n" ] } ], "prompt_number": 18 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Verify that the above vectors are indeed stationary distributions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hint: Use [np.dot](http://docs.scipy.org/doc/numpy/reference/generated/numpy.dot.html)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Write your own code\n", "svec1=(1,0,0,0,0,0)\n", "svec2=(0,1/3,0,0,2/3,0)\n", "np.dot(svec1,P)\n", "np.dot(svec2,P)\n", "print (np.dot(svec1,P))\n", "print (np.dot(svec2,P))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 1. 0. 0. 0. 0. 0.]\n", "[ 0. 0.33333333 0. 0. 0.66666667 0. ]\n" ] } ], "prompt_number": 43 }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Simulation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us simulate our Markov chain `mc0`.\n", "The `simualte` method generates a sample path from an initial state as specified by the `init` argument,\n", "of length specified by `sample_size`, which is set to 1000 by default when omitted." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A sample path from state `0`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "mc0.simulate(init=0, sample_size=50)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 6, "text": [ "array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n", " 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n", " 0, 0, 0, 0])" ] } ], "prompt_number": 6 }, { "cell_type": "markdown", "metadata": {}, "source": [ "As is clear from the transition matrix `P`,\n", "if it starts at state `0`, the chain stays there forever,\n", "i.e., `0` is an absorbing state, a state that constitutes a singleton recurrent class." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Start with state `1`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "mc0.simulate(init=1, sample_size=50)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 7, "text": [ "array([1, 4, 1, 4, 1, 4, 1, 4, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 4, 1, 4, 1,\n", " 4, 4, 4, 1, 4, 1, 4, 1, 4, 4, 1, 4, 1, 4, 1, 4, 1, 4, 4, 4, 4, 1, 4,\n", " 1, 4, 1, 4])" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can observe that the chain stays in the recurrent class $\\{1, 4\\}$ and\n", "visits states `1` and `4` with certain frequencies." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us compute the frequency distribution along a sample path.\n", "We will repeat this operation, so let us write a function for it." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def time_series_dist(mc, init, ts_length=100):\n", " \"\"\"\n", " Return a distribution of visits by a sample path of length ts_length\n", " of mc with an initial state init.\n", " \n", " \"\"\"\n", " X = mc.simulate(init=init, sample_size=ts_length)\n", " bins = np.arange(mc.n+1)\n", " hist, bin_edges = np.histogram(X, bins=bins)\n", " dist = hist/len(X)\n", " return dist" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 8 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is a frequency distribution along a sample path from initial state `1`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "time_series_dist(mc0, init=1)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 9, "text": [ "array([ 0. , 0.35, 0. , 0. , 0.65, 0. ])" ] } ], "prompt_number": 9 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us visualize the distribution." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def draw_histogram(distribution,\n", " ax=None, title=None, xlabel=None, ylabel=None, ylim=(0, 1)):\n", " \"\"\"\n", " Plot the given distribution.\n", " \n", " \"\"\"\n", " if ax is None:\n", " fig, ax = plt.subplots()\n", " n = len(distribution)\n", " ax.bar(np.arange(n), distribution, align='center')\n", " ax.set_xlim(-0.5, (n-1)+0.5)\n", " ax.set_ylim(*ylim)\n", " if title:\n", " ax.set_title(title)\n", " if xlabel:\n", " ax.set_xlabel(xlabel)\n", " if ylabel:\n", " ax.set_ylabel(ylabel)\n", " if ax is None:\n", " plt.show()" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 11 }, { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 11 }, { "cell_type": "code", "collapsed": false, "input": [ "init = 1\n", "draw_histogram(time_series_dist(mc0, init=init),\n", " title='Time series distribution with init={0}'.format(init),\n", " xlabel='States')\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAAEZCAYAAABYR6TIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFx1JREFUeJzt3XuYZHV95/H3xwG8IKMoLoaL4mW8kPUWXcQ1aBuMjkZl\nkycuDkGXyK5kdzFZjQmLyWpn2ayarNEk+iC6BC8xjIm6LhoErw2IiuIChjBjIIgyDI4IBBWIzjjf\n/eOcxqKn+jZUT/Wv+/16nn7mnDq/+p3vqer51K9+51R1qgpJUlvuNe4CJEmLZ3hLUoMMb0lqkOEt\nSQ0yvCWpQYa3JDXI8B6zJFcmeda461iMJOcmefkS72MyyQf65Ycl+UGSjKjv05P8fr88keT6UfTb\n93dUks2j6m/U+09yWJKdSRb0f3/wsRplW91z8TrvpZXkh8D0g7wv8M/AT/r1V1XV2WMpbJlL8kbg\n0VW14BeJJCcAJ1bVUYu4zwTwgao6dNFFdvffSVfntbtz/6WW5DrglVX1uX79MOBaYK+q2rmE+51g\nNx/XJP8SeCvwc8CDq8pB5hA+KEusqu5fVftV1X7At4AXTa+3FtzpjbuOe2KhI87FdrsEfY5Ksbzr\nG+bHwEbgxHEXspwZ3mOW5Lokv9AvTyb5myQfSPL9JF9Psi7JqUm2JflWkl8cuO8DkpyZZGuSLUlO\nmy2ckhyR5NIktyX5TpK3Dmw7MskXk9ya5PIkzx7YNpXkfyS5GPgh8Mj+thMH2rwyyVVJbklyXpKH\nDWx7W1/7bf3x/Ows9T0iyQX9cX8KOGBg293e6ic5Ick/9m2vTXJckscB7wKe0U+x3NK3fW//dv7c\n/l3Qc/rbTpux/1OT3JTkm0mOm3H8g8d6QpKL+uUL+5uv6Pf50pnTMEke3/dxaz9F9uKBbe9N8s4k\nn+iP5ctJHjnL4/O+JK/tlw/uH4//1K8/KsnN/fJd+0837fQw4ON9fa8b6PL4/vfppiSvH7bPgRpP\nG+h7S5LX9s/p1nTvdu7WNsn9gE8CB/X7/X6Sh862j5mq6h+q6izgqoXeZzUyvMdv5rzVi4D3A/sD\nlwGf7m8/CDgNOGOg7XvpRimPAp4CPA/497Ps50+Bt1XVA4BHAn8NXRAAnwD+e1XtD7wO+EiSBw/c\n9/i+3+l3DzVdd5JjgFOBX6YL3IuAs/ttzweOAtb1+30pcPMs9f0V8FXgwf1x/rshjw1J9u2PZX1V\nrQWeAVxeVZuBk4Av9e9qHjRwtw3AaVV1f+ALg/X3Htrv96B+v+9Osq7fNrPtXapq+lzFE/t9/s2M\nWvcGPg6cBzwEeDXwwSSPGWh2LDBJ93xfA/zh0EcHpoCJfvnZdFMfzxpYv3DmHfopp2/z03d7/2tg\n8zOBxwBHA2/oX/yGHiZ3P/4DgbV0j9WJwDuTPGCwbVXdAawHtvb7XVtV3+lfZG+d5eeWJIfMUoOG\nMLyXnwur6tNV9RPgw3Sh8uZ+/UPAYUnWJjkQeAHwmqq6s6puAt4OvGyWfn8MrEtyQFXdUVWX9Lcf\nD5xbVecBVNVngEuBX+q3F/DeqtpUVTuraseMfn8DeFNVfaOfQ30T8OR+9P1jusB/fJJ79W2+M7Ow\nvu3TgP9WVdur6iK60Jvt7f5O4AlJ7ltV26pqeoQ2rH0BH6uqL/XH96NZ2k7v+0Lgb+lC9Z46Eti3\nqt5cVTuq6vN0L5QbBtp8tKou7Z/fDwJPnqWvC4GfTxK6F8Q/ogtg6ML7gkXW9gdV9aOq+jpwBfCk\nOdoOPlbb6V7of1JVn6R7N/bYIW13eS6q6q+qav9Zfh5UVVsWeQyrmuG9/Hx3YPlO4Hv107PKd/b/\n3h94OLA3cOP06IVu2uAhs/R7It1Ia1OSrySZDueHAy8dHAXRhcLg29y5rsZ4OPCnA/edHlkf1IfV\nO4B3AtuSnJFkvyF9HATcWlV3Dtz2rWE7q6rb6YL1N4Ct/ZTDY4e1XWD9zLLvn5nnPgtx0JB9f6u/\nHboXlm0D2+6ke253UVX/CNxOF+5H0b0IbO1H8c9i8eE9+CJ6B93J9IW4ecaJzjtmq1lLy/Bu1/XA\nj+jOxk+PXh5QVU8Y1riqrqmq46rqIcBbgA/3c5PfprsqYHAUtF9V/dHg3eeo49t0V80M3n/fqvpy\nv98/r6qnAYfTvXj8zpA+bgT27+uZ9vDZ9ltVn6qq59G9wGwG3rOAOnfpZmB52L639su3c/dgW/Dc\nbd/Hof1oebDvGxbRx6AL6Kae9q6qrf36CXRTLpfPcp9RXE62O4/rsCmvX+vnwIf9fN9pk8UxvBtV\nVTcCnwL+JMl+Se7Vn7gaes14kuOTTI/Kb6P7z/UT4C+BFyd5XpI1Se7Tn5g6ePDuc5TyLuD1SQ7v\n9/OAJC/tl5+W5On93O8d3P0yycFj+RbdVM0fJNk7yc/Tzf0PO45/keSYfu57O124Tve5DTik399c\ntWfI7dP7Popuymh6/vpy4FeS3DfJo9n1CohtdOcchrmE7rh/t+97oj+ujXPUNpcLgJP56fz2VL9+\n0cC7s5nmqm/QbLUMe6zm6mO67TbgwUnWTm+sqg8OXGk182ft4LRJkvsA+/TL905y7wXWsGoY3svL\nsJNjc62/gu4X/CrgFrrAmW1k+HzgyiQ/AN4GvKyf89wCHAO8nm7K5tvAb3P3/7Czjryq6mN0I/mN\nSW4D/q7fF3Qntt7d13Yd8D3gj2fp6jjg6X3bNwDvm+W47wW8hm70ejPdFMJ/7Ld9Fvh74DtJvjtw\nv2GP4eBtNwK30o2UPwCcVFX/0G97G93c/TbgLLoXu8H7TgLv66eNfnWw76r6MfBiunMTN9FNIb18\noO+FPN+DLqSbopgO74uB+7LrycrBPt4E/H5f32vn2Mds+51Z41z1DR77ZroT19f2JyMX/I4l3bXo\ndwBX9v3dCWxa6P1Xi3k/pJPkL+hGIt+d7S15kj+j+wW9Azihqi4bdaGSpJ9ayMj7LLrLfoZK8kK6\nT5itA14FnD6i2iRJs5g3vPvLtm6do8lL6N/i9pefPbC/jE2StERGMed9MHe/HGoL4FljSVpCozph\nOfNstN92JUlLaK8R9HEDMPjNYYcw5DrWJAa6JO2Gqtrlcs1RjLzPobtkjSRHAv9UVduGNayqPfbz\nxje+cY/ub0//rOTjW8nH5vG1/7Onj2828468k5xN990JB6T7trI30n0sm6o6o6rOTfLCJNfQfWDi\n1xcZ/pKkRZo3vKtqwwLanDyaciRJC7FiP2E5MTEx7hKW1Eo+vpV8bODxtW65HN8e+zNoSWpP7UuS\nVook1BKdsJQk7WGGtyQ1yPCWpAYZ3pLUIMNbkhpkeEtSgwxvSWqQ4S1JDTK8JalBhrckNcjwlqQG\nGd6S1CDDW5IaZHhLUoMMb0lqkOEtSQ0yvCWpQYa3JDXI8JakBhnektQgw1uSGmR4S1KDDG9JapDh\nLUkNMrwlqUGGtyQ1yPCWpAYZ3pLUIMNbkhpkeEtSgwxvSWqQ4S1JDTK8JalBhrckNWje8E6yPsnm\nJFcnOWXI9gOSnJfk8iRXJjlhSSqVJN0lVTX7xmQN8A3gucANwFeBDVW1aaDNJHDvqjo1yQF9+wOr\naseMvmqufUmSdpWEqsrM2+cbeR8BXFNV11XVdmAjcMyMNjcCa/vltcDNM4NbkjRae82z/WDg+oH1\nLcDTZ7R5D/C5JFuB/YB/O7ryJEnDzBfeC5nneD1weVVNJHkU8OkkT6qqH8xsODk5edfyxMQEExMT\niyhVkla+qakppqam5m0335z3kcBkVa3v108FdlbVWwbanAv8YVVd3K9/Fjilqi6d0Zdz3pK0SLs7\n530psC7JYUn2AY4FzpnRZjPdCU2SHAg8Frj2npcsSZrNnNMmVbUjycnA+cAa4Myq2pTkpH77GcD/\nBM5KcgXdi8HvVtUtS1y3JK1qc06bjHRHTptI0qLt7rSJJGkZMrwlqUGGtyQ1yPCWpAYZ3pLUIMNb\nkhpkeEtSgwxvSWqQ4S1JDTK8JalBhrckNcjwlqQGGd6S1CDDW5IaZHhLUoMMb0lqkOEtSQ0yvCWp\nQYa3JDXI8JakBhnektQgw1uSGmR4S1KDDG9JapDhLUkNMrwlqUGGtyQ1yPCWpAYZ3pLUoL3GXYCk\n8Ugy7hIWrarGXcKyYXhLq1pLYdjei81SctpEkhpkeEtSgwxvSWqQ4S1JDZo3vJOsT7I5ydVJTpml\nzUSSy5JcmWRq5FVKku4mc116k2QN8A3gucANwFeBDVW1aaDNA4GLgedX1ZYkB1TV94b0VV7mIy0f\n3aWCLf2fzKq8VDAJVbXLpTbzjbyPAK6pquuqajuwEThmRpvjgI9U1RaAYcEtSRqt+cL7YOD6gfUt\n/W2D1gEPSvL5JJcmefkoC5Qk7Wq+D+ks5D3K3sDPAUcD9wO+lOTLVXX1PS1OkjTcfOF9A3DowPqh\ndKPvQdcD36uqO4E7k1wIPAnYJbwnJyfvWp6YmGBiYmLxFUvSCjY1NcXU1NS87eY7YbkX3QnLo4Gt\nwFfY9YTl44B3AM8H7g1cAhxbVVfN6MsTltIy4gnLNsx2wnLOkXdV7UhyMnA+sAY4s6o2JTmp335G\nVW1Och7wdWAn8J6ZwS1JGq05R94j3ZEjb2lZceTdht29VFCStAwZ3pLUIMNbkhpkeEtSgwxvSWqQ\n4S1JDTK8JalBhrckNcjwlqQGGd6S1CDDW5IaZHhLUoMMb0lqkOEtSQ0yvCWpQYa3JDXI8JakBhne\nktQgw1uSGmR4S1KDDG9JapDhLUkNMrwlqUGGtyQ1yPCWpAYZ3pLUIMNbkhpkeEtSgwxvSWqQ4S1J\nDTK8JalBhrckNcjwlqQGGd6S1CDDW5IaZHhLUoPmDe8k65NsTnJ1klPmaPevkuxI8iujLVGSNNOc\n4Z1kDfAOYD1wOLAhyeNnafcW4DwgS1CnJGnAfCPvI4Brquq6qtoObASOGdLu1cCHgZtGXJ8kaYj5\nwvtg4PqB9S39bXdJcjBdoJ/e31Qjq06SNNR84b2QIH478F+rquimTJw2kaQlttc8228ADh1YP5Ru\n9D3oqcDGJAAHAC9Isr2qzpnZ2eTk5F3LExMTTExMLL5iSVrBpqammJqamrddugHzLBuTvYBvAEcD\nW4GvABuqatMs7c8CPl5VHx2yrebal6Q9qxtwtfR/MqzGDElCVe0yozHnyLuqdiQ5GTgfWAOcWVWb\nkpzUbz9jSaqVJM1pzpH3SHfkyFtaVhx5t2G2kbefsJSkBhnektQgw1uSGmR4S1KDDG9JapDhLUkN\nMrwlqUHzfTxe90D/lQHNWY3X0kqtMbyXXGtB2OYLjrTaOG0iSQ0yvCWpQYa3JDXI8JakBhnektQg\nw1uSGmR4S1KDDG9JapDhLUkNMrwlqUGGtyQ1yPCWpAYZ3pLUIMNbkhpkeEtSgwxvSWqQ4S1JDTK8\nJalBhrckNcjwlqQGGd6S1CDDW5IaZHhLUoMMb0lqkOEtSQ0yvCWpQYa3JDVoQeGdZH2SzUmuTnLK\nkO2/luSKJF9PcnGSJ46+VEnStHnDO8ka4B3AeuBwYEOSx89odi3wrKp6InAa8O5RFypJ+qmFjLyP\nAK6pquuqajuwEThmsEFVfamqbutXLwEOGW2ZkqRBCwnvg4HrB9a39LfN5kTg3HtSlCRpbnstoE0t\ntLMkzwFeCTxz2PbJycm7licmJpiYmFho15K0KkxNTTE1NTVvu1TNnc1JjgQmq2p9v34qsLOq3jKj\n3ROBjwLrq+qaIf3UfPtaaZKwiNe+ZSKstudptWrv93N1/m4moaoy8/aFTJtcCqxLcliSfYBjgXNm\ndP4wuuA+flhwS5JGa95pk6rakeRk4HxgDXBmVW1KclK//QzgDcD+wOndqznbq+qIpStbkla3eadN\nRrYjp00asTrfmq5G7f1+rs7fzXsybSJJWmYMb0lqkOEtSQ0yvCWpQYa3JDXI8JakBhnektQgw1uS\nGmR4S1KDDG9JapDhLUkNMrwlqUGGtyQ1yPCWpAYZ3pLUIMNbkhpkeEtSgwxvSWqQ4S1JDTK8JalB\nhrckNcjwlqQGGd6S1CDDW5IaZHhLUoMMb0lqkOEtSQ0yvCWpQYa3JDXI8JakBhnektQgw1uSGmR4\nS1KDDG9JapDhLUkNmje8k6xPsjnJ1UlOmaXNn/Xbr0jylNGXKUkaNGd4J1kDvANYDxwObEjy+Blt\nXgg8uqrWAa8CTl+iWhdlampq3CUssalxF7BkVvpzt9KPbyX/bsLyef7mG3kfAVxTVddV1XZgI3DM\njDYvAd4HUFWXAA9McuDIK12k5fIAL52pcRewZFb6c7fSj28l/27C8nn+5gvvg4HrB9a39LfN1+aQ\ne16aJGk284V3LbCf7Ob9JEm7IVWz52ySI4HJqlrfr58K7Kyqtwy0eRcwVVUb+/XNwLOratuMvgx0\nSdoNVTVzgMxe89znUmBdksOArcCxwIYZbc4BTgY29mH/TzODe7adS5J2z5zhXVU7kpwMnA+sAc6s\nqk1JTuq3n1FV5yZ5YZJrgNuBX1/yqiVplZtz2kSStDytuE9YLuRDRa1K8hdJtiX5u3HXshSSHJrk\n80n+PsmVSX5z3DWNUpL7JLkkyeVJrkrypnHXNGpJ1iS5LMnHx13LUkhyXZKv98f4lbHWspJG3v2H\nir4BPBe4AfgqsKGqNo21sBFJchTwQ+D9VfWEcdczakkeCjy0qi5Pcn/ga8C/WSnPH0CS+1XVHUn2\nAr4AvK6qvjDuukYlyWuBpwL7VdVLxl3PqCX5JvDUqrpl3LWstJH3Qj5U1Kyqugi4ddx1LJWq+k5V\nXd4v/xDYBBw03qpGq6ru6Bf3oTuPNPYQGJUkhwAvBP43u14+vJIsi2NbaeG9kA8VqQH9FU5PAS4Z\nbyWjleReSS4HtgGfr6qrxl3TCL0N+B1g57gLWUIFfCbJpUn+wzgLWWnhvXLmgFaxfsrkw8Bv9SPw\nFaOqdlbVk+k+hfysJBNjLmkkkrwI+G5VXcYyGZkukWdW1VOAFwD/uZ/KHIuVFt43AIcOrB9KN/pW\nI5LsDXwE+Muq+ti461kqVXUb8LfA08Zdy4j8a+Al/Zzw2cAvJHn/mGsauaq6sf/3JuD/0E3VjsVK\nC++7PlSUZB+6DxWdM+aatEBJApwJXFVVbx93PaOW5IAkD+yX7wv8InDZeKsajap6fVUdWlWPAF4G\nfK6qXjHuukYpyf2S7Ncv7ws8DxjblV8rKryragfdpz3PB64CPrTCrlQ4G/gi8Jgk1ydZaR+IeiZw\nPPCc/lKsy5KsH3dRI/QzwOf6Oe9LgI9X1WfHXNNSWYlTmAcCFw08f5+oqk+Nq5gVdamgJK0WK2rk\nLUmrheEtSQ0yvCWpQYa3JDXI8JakBhnektQgw1srRpLf679K9or+GvEjkvxW/4GY+e77XxbSTlou\nvM5bK0KSZwBvpfv7qduTPAi4D3Ax8LSqunme+39zIe2k5cKRt1aKhwLf678KmP77ln+V7itlP5/k\nswBJTk/y1X6EPtnf9ptD2j0vyReTfC3JX/cfhybJm/s/FnFFkj/e40cp9Rx5a0Xow/ULwP2Az9B9\nNcKFM788P8n+VXVr/4c7PgO8uqquHGyX5AC6L8daX1V39n+RaR/gncAXq+pxfV9rq+r7e/xgJRx5\na4Woqtvp/oLLq4CbgA8lOaHfPPgVpccm+Rrw/4CfBQ4f0t2R/e1fTHIZ8ArgYcBtwD8nOTPJLwN3\nLsWxSAsx51+Pl1pSVTuBC4AL+r/zecL0JoAkjwB+m25u+7YkZ9HNiw/z6ao6buaNSY4Ajqabkjm5\nX5b2OEfeWhGSPCbJuoGbngJcB/wAWNvftha4Hfh+kgPpvlB/2mC7S4BnJnlU3/e+Sdb1UzMPrKpP\nAq8FnrRUxyPNx5G3Vor7A3/ef1/2DuBquimU44DzktxQVUf30yCb6f5c3uAf/n33jHYnAGcnuXe/\n/ffoAv7/JrkP3VTMa/bEgUnDeMJSkhrktIkkNcjwlqQGGd6S1CDDW5IaZHhLUoMMb0lqkOEtSQ0y\nvCWpQf8fBho4mtLqCvwAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 12 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Observe that the distribution is close to the stationary distribution\n", "`[0, 1/3, 0, 0, 2/3, 0]`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Start with state `2`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "init = 2\n", "draw_histogram(time_series_dist(mc0, init=init),\n", " title='Time series distribution with init={0}'.format(init),\n", " xlabel='States')\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAAEZCAYAAABYR6TIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF3xJREFUeJzt3XuYZHV95/H3xwG8IKMoLoaL4mW8kMdbdBGXoG0wOroq\niU8MDkFDdFeyuxhXc2EhrrTLJmouq0n0QTTEWxRM1HVREbw2ICqKCxhkxoUgyjA4IhBQIDrjfPeP\ncxqKnurLMNVT/et+v55nnjmnzq9+53uquj/1q985VZ2qQpLUlnuNuwBJ0o4zvCWpQYa3JDXI8Jak\nBhnektQgw1uSGmR4j1mSy5M8c9x17IgkZyd5+SLvYzLJB/vlhyX5cZKMqO9Tk7yhX55Icu0o+u37\nOzzJhlH1N+r9JzkoybYkC/rdH3ysRtlWOy9e5724kvwEmH6Q9wT+Ffh5v/7qqjpjLIUtcUlOBh5d\nVQt+kUhyLPCqqjp8B+4zAXywqg7c4SK7+2+jq/Pqe3L/xZbkGuCVVfXFfv0g4Gpgt6ratoj7neAe\nPq5Jfht4DbAGuBX4MHBSVf18zjuuMI68F1lV3b+q9qqqvYDvAS+cXm8tuNMbdx07Y6Ejzh3tdhH6\nHJViadc3zH2B1wIPBp4OHAH8wVgrWoIM7zFLck2SX+mXJ5P8Y5IPJrk1ybeSrElyYpLNSb6X5FcH\n7vuAJKcn2ZRkY5JTZgunJIckuTjJLUl+kOQvB7YdmuQrSW5OcmmSZw1sm0ryP5NcCPwEeGR/26sG\n2rwyyRVJbkpyTpKHDWx7W1/7Lf3x/OIs9T0iyXn9cX8W2Gdg293e6ic5Nsk/922vTnJ0kscB7wKe\n0U+x3NS3fV//dv7s/l3Qs/vbTpmx/xOT3JDku0mOnnH8g8d6bJIL+uXz+5sv6/f50pnTMEke3/dx\ncz9F9qKBbe9L8s4kn+qP5WtJHjnL4/P+JK/vl/fvH4//3K8/KsmN/fKd+0837fQw4JN9fYMBeEz/\n83RDkpOG7XOgxlMG+t6Y5PX9c7op3budu7VNcj/gM8B+/X5vTfLQ2fYxU1W9q6ourKqtVbUJ+BBw\n2ELvv1IY3uM3c97qhcAHgL2BS4DP9bfvB5wCnDbQ9n3Az4BHAU8Bngv8h1n281fA26rqAcAjgX+A\nLgiATwH/o6r2phvhfCzJgwfue0zf7/S7h5quO8mRwInAr9MF7gXAGf225wGHA2v6/b4UuHGW+j4M\nfINutHUK8NtDHhuS7Nkfy9qqWg08A7i0qjYAxwFf7d/VPGjgbuuAU6rq/sCXB+vvPbTf7379ft+d\nZE2/bWbbO1XV9LmKJ/b7/McZte4OfBI4B3gI3VTAh5I8ZqDZUcAk3fN9FfAnQx8dmAIm+uVn0U19\nPHNg/fyZd+innL7PXe/2/mJg82HAY+hGtW/sX/yGHiZ3P/59gdV0j9WrgHcmecBg26q6HVgLbOr3\nu7qqftC/yN48y7+bkhwwSw3PAi6fZduKZXgvPedX1ef6+b2P0oXKW/r1jwAHJVmdZF/g+cDrquqO\nqroBeDvwsln6/RmwJsk+VXV7VV3U334McHZVnQNQVZ8HLgb+fb+9gPdV1fqq2lZVW2f0+7vAm6vq\nO/0c6puBJ/ej75/RBf7jk9yrb/ODmYX1bZ8G/Peq2lJVF9CF3mxv97cBT0hy36raXFVXTHc1pG0B\nn6iqr/bH99NZ2k7v+3zg03ShurMOBfasqrf0o8gv0b1Qrhto8/Gqurh/fj8EPHmWvs4HfjlJ6F4Q\n/4y7RqPPAs7bwdreVFU/rapvAZcBT5qj7eBjtYXuhf7nVfUZundjjx3Sdrvnoqo+XFV7z/LvQVW1\ncbsdJ68Efgn4i5nbVjrDe+n54cDyHcCP6q6zynf0/98feDiwO3D99OiFbtrgIbP0+yq6kdb6JF9P\nMh3ODwdeOjgKoguFwbe5c12N8XDgrwbuOz2y3q8Pq3cA7wQ2JzktyV5D+tgPuLmq7hi47XvDdlZV\nt9EF6+8Cm/oph8cOa7vA+pll378wz30WYr8h+/5efzt0LyybB7bdQffcbqeq/hm4jS7cD6d7EdjU\nj+KfyY6H9+CL6O10J9MX4sYZJzpvn63mnZXk14A/BZ5fVTctxj5aZni361rgp8CDB0YvD6iqJwxr\nXFVXVdXRVfUQ4K3AR/u5ye/TXRUwOAraq6r+bPDuc9TxfbqrZgbvv2dVfa3f799U1dOAg+lePP5w\nSB/XA3v39Ux7+Gz7rarPVtVz6V5gNgDvWUCd23UzsDxs35v65du4e7AteO627+PAfrQ82Pd1O9DH\noPPopp527+eCzwOOpZtyuXSW+4zicrJ78rgOm/L6rX4OfNi/WwenTZKsBd5NN+Xz7Z0rf3kyvBtV\nVdcDnwX+V5K9ktyrP3E19JrxJMckmR6V30L3y/Vz4O+BFyV5bpJVSe7Tn5jaf/Duc5TyLuCkJAf3\n+3lAkpf2y09L8vR+7vd27n6Z5OCxfI9uquZNSXZP8st0c//DjuPfJDmyn/veQheu031uBg7o9zdX\n7Rly+/S+D6ebMpqev74UeEmS+yZ5NN07mEGb6c45DHMR3XH/Ud/3RH9cZ85R21zOA47nrvntqX79\ngoF3ZzPNVd+g2WoZ9ljN1cd0283Ag5Osnt5YVR8auNJq5r/V09Mm6U7gfwh4SVVdvMB9rziG99Iy\n7OTYXOuvAPYArgBuoguc2UaGzwMuT/Jj4G3Ay/o5z43AkcBJdFM23wd+n7v/ws468qqqT9CN5M9M\ncgvwT/2+oDux9e6+tmuAHwF/PktXR9NdFnYT8Ebg/bMc972A19GNXm+km0L4T/22LwDfBn6Q5IcD\n9xv2GA7edj1wM91I+YPAcVX1//ptb6Obu98MvJfuxW7wvpPA+/tpo98Y7Luqfga8iO7cxA10U0gv\nH+h7Ic/3oPPppiimw/tCusvqZp6sHOzjzcAb+vpeP8c+ZtvvzBrnqm/w2DfQnbi+uj8ZuSPvWN5A\nd67kMwMj80/vwP1XhHk/pJPk7+hGIj+c7S15kr+m+wG9HTi2qi4ZdaGSpLssZOT9XrrLfoZK8gK6\nT5itAV4NnDqi2iRJs5g3vPvLtm6eo8mL6d/i9pefPbC/jE2StEhGMee9P3e/HGojMNvF9pKkERjV\nCcuZZ6P9titJWkS7jaCP64DBbw47gCHXsSYx0CXpHqiq7S7XHMXI+yy6S9ZIcijwL1W1eVjDqtpl\n/04++eRdur9d/W85H99yPjaPr/1/u/r4ZjPvyDvJGXTfnbBPum8rO5nuY9lU1WlVdXaSFyS5iu4D\nE7+zg+EvSdpB84Z3Va1bQJvjR1OOJGkhlu0nLCcmJsZdwqJazse3nI8NPL7WLZXj22V/Bi1J7ap9\nSdJykYRapBOWkqRdzPCWpAYZ3pLUIMNbkhpkeEtSgwxvSWqQ4S1JDTK8JalBhrckNcjwlqQGGd6S\n1CDDW5IaZHhLUoMMb0lqkOEtSQ0yvCWpQYa3JDXI8JakBhnektQgw1uSGmR4S1KDDG9JapDhLUkN\nMrwlqUGGtyQ1yPCWpAYZ3pLUIMNbkhpkeEtSgwxvSWqQ4S1JDTK8JalBhrckNcjwlqQGzRveSdYm\n2ZDkyiQnDNm+T5Jzklya5PIkxy5KpZKkO6WqZt+YrAK+AzwHuA74BrCuqtYPtJkE7l1VJybZp2+/\nb1VtndFXzbUvSdL2klBVmXn7fCPvQ4CrquqaqtoCnAkcOaPN9cDqfnk1cOPM4JYkjdZu82zfH7h2\nYH0j8PQZbd4DfDHJJmAv4DdHV54kaZj5wnsh8xwnAZdW1USSRwGfS/KkqvrxzIaTk5N3Lk9MTDAx\nMbEDpUrS8jc1NcXU1NS87eab8z4UmKyqtf36icC2qnrrQJuzgT+pqgv79S8AJ1TVxTP6cs5bknbQ\nPZ3zvhhYk+SgJHsARwFnzWizge6EJkn2BR4LXL3zJUuSZjPntElVbU1yPHAusAo4varWJzmu334a\n8KfAe5NcRvdi8EdVddMi1y1JK9qc0yYj3ZHTJpK0w+7ptIkkaQkyvCWpQYa3JDXI8JakBhnektQg\nw1uSGmR4S1KDDG9JapDhLUkNMrwlqUGGtyQ1yPCWpAYZ3pLUIMNbkhpkeEtSgwxvSWqQ4S1JDTK8\nJalBhrckNcjwlqQGGd6S1CDDW5IaZHhLUoMMb0lqkOEtSQ0yvCWpQYa3JDXI8JakBhnektQgw1uS\nGmR4S1KDdht3AZLGI8m4S9hhVTXuEpYMw1ta0VoKw/ZebBaT0yaS1CDDW5IaNG94J1mbZEOSK5Oc\nMEubiSSXJLk8ydTIq5Qk3U3mOgGQZBXwHeA5wHXAN4B1VbV+oM0DgQuB51XVxiT7VNWPhvRVnmyQ\nlo7uhGVLv5NZkScsk1BV2034zzfyPgS4qqquqaotwJnAkTPaHA18rKo2AgwLbknSaM0X3vsD1w6s\nb+xvG7QGeFCSLyW5OMnLR1mgJGl7810quJD3KLsDvwQcAdwP+GqSr1XVlTtbnCRpuPnC+zrgwIH1\nA+lG34OuBX5UVXcAdyQ5H3gSsF14T05O3rk8MTHBxMTEjlcsScvY1NQUU1NT87ab74TlbnQnLI8A\nNgFfZ/sTlo8D3gE8D7g3cBFwVFVdMaMvT1hKS4gnLNsw2wnLOUfeVbU1yfHAucAq4PSqWp/kuH77\naVW1Ick5wLeAbcB7Zga3JGm05hx5j3RHjrylJcWRdxvu6aWCkqQlyPCWpAYZ3pLUIMNbkhpkeEtS\ngwxvSWqQ4S1JDTK8JalBhrckNcjwlqQGGd6S1CDDW5IaZHhLUoMMb0lqkOEtSQ0yvCWpQYa3JDXI\n8JakBhnektQgw1uSGmR4S1KDDG9JapDhLUkNMrwlqUGGtyQ1yPCWpAYZ3pLUIMNbkhpkeEtSgwxv\nSWqQ4S1JDTK8JalBhrckNcjwlqQGGd6S1CDDW5IaNG94J1mbZEOSK5OcMEe7f5tka5KXjLZESdJM\nc4Z3klXAO4C1wMHAuiSPn6XdW4FzgCxCnZKkAfONvA8Brqqqa6pqC3AmcOSQdq8BPgrcMOL6JElD\nzBfe+wPXDqxv7G+7U5L96QL91P6mGll1kqSh5gvvhQTx24H/VlVFN2XitIkkLbLd5tl+HXDgwPqB\ndKPvQU8FzkwCsA/w/CRbquqsmZ1NTk7euTwxMcHExMSOVyxJy9jU1BRTU1Pztks3YJ5lY7Ib8B3g\nCGAT8HVgXVWtn6X9e4FPVtXHh2yrufYladfqBlwt/U6GlZghSaiq7WY05hx5V9XWJMcD5wKrgNOr\nan2S4/rtpy1KtZKkOc058h7pjhx5S0uKI+82zDby9hOWktQgw1uSGmR4S1KDDG9JapDhLUkNMrwl\nqUGGtyQ1aL6Px2sn9F8Z0JyVeC2t1BrDe9G1FoRtvuBIK43TJpLUIMNbkhpkeEtSgwxvSWqQ4S1J\nDTK8JalBhrckNcjwlqQGGd6S1CDDW5IaZHhLUoMMb0lqkOEtSQ0yvCWpQYa3JDXI8JakBhnektQg\nw1uSGmR4S1KDDG9JapDhLUkNMrwlqUGGtyQ1yPCWpAYZ3pLUIMNbkhpkeEtSgxYU3knWJtmQ5Mok\nJwzZ/ltJLkvyrSQXJnni6EuVJE2bN7yTrALeAawFDgbWJXn8jGZXA8+sqicCpwDvHnWhkqS7LGTk\nfQhwVVVdU1VbgDOBIwcbVNVXq+qWfvUi4IDRlilJGrSQ8N4fuHZgfWN/22xeBZy9M0VJkua22wLa\n1EI7S/Js4JXAYcO2T05O3rk8MTHBxMTEQruWpBVhamqKqampedulau5sTnIoMFlVa/v1E4FtVfXW\nGe2eCHwcWFtVVw3pp+bb13KThB147Vsiwkp7nlaq9n4+V+bPZhKqKjNvX8i0ycXAmiQHJdkDOAo4\na0bnD6ML7mOGBbckabTmnTapqq1JjgfOBVYBp1fV+iTH9dtPA94I7A2c2r2as6WqDlm8siVpZZt3\n2mRkO3LapBEr863pStTez+fK/NncmWkTSdISY3hLUoMMb0lqkOEtSQ0yvCWpQYa3JDXI8JakBhne\nktQgw1uSGmR4S1KDDG9JapDhLUkNMrwlqUGGtyQ1yPCWpAYZ3pLUIMNbkhpkeEtSgwxvSWqQ4S1J\nDTK8JalBhrckNcjwlqQGGd6S1CDDW5IaZHhLUoMMb0lqkOEtSQ0yvCWpQYa3JDXI8JakBhnektQg\nw1uSGmR4S1KDDG9JatC84Z1kbZINSa5McsIsbf66335ZkqeMvkxJ0qA5wzvJKuAdwFrgYGBdksfP\naPMC4NFVtQZ4NXDqItW6Q6ampsZdwiKbGncBi2a5P3fL/fiW888mLJ3nb76R9yHAVVV1TVVtAc4E\njpzR5sXA+wGq6iLggUn2HXmlO2ipPMCLZ2rcBSya5f7cLffjW84/m7B0nr/d5tm+P3DtwPpG4OkL\naHMAsHmnq9OSlmTR+n7Tm960aH1X1aL1Le0q8428F/pTPvO32N+OFaMW4d/Ji9SvP5ZaPjLXKCTJ\nocBkVa3t108EtlXVWwfavAuYqqoz+/UNwLOqavOMvvzNkaR7oKq2e5s737TJxcCaJAcBm4CjgHUz\n2pwFHA+c2Yf9v8wM7tl2Lkm6Z+YM76ramuR44FxgFXB6Va1Pcly//bSqOjvJC5JcBdwG/M6iVy1J\nK9yc0yaSpKVp2X3CciEfKmpVkr9LsjnJP427lsWQ5MAkX0ry7SSXJ/m9cdc0Sknuk+SiJJcmuSLJ\nm8dd06glWZXkkiSfHHctiyHJNUm+1R/j18day3IaefcfKvoO8BzgOuAbwLqqWj/WwkYkyeHAT4AP\nVNUTxl3PqCV5KPDQqro0yf2BbwK/tlyeP4Ak96uq25PsBnwZ+IOq+vK46xqVJK8HngrsVVUvHnc9\no5bku8BTq+qmcdey3EbeC/lQUbOq6gLg5nHXsViq6gdVdWm//BNgPbDfeKsaraq6vV/cg+480thD\nYFSSHAC8APhbtr98eDlZEse23MJ72AeG9h9TLdoJ/RVOTwEuGm8lo5XkXkkupfsQ25eq6opx1zRC\nbwP+ENg27kIWUQGfT3Jxkv84zkKWW3gvnzmgFayfMvko8Np+BL5sVNW2qnoy3aeQn5lkYswljUSS\nFwI/rKpLWCIj00VyWFU9BXg+8F/6qcyxWG7hfR1w4MD6gXSjbzUiye7Ax4C/r6pPjLuexVJVtwCf\nBp427lpG5N8BL+7nhM8AfiXJB8Zc08hV1fX9/zcA/5tuqnYsllt43/mhoiR70H2o6Kwx16QFSvdl\nKacDV1TV28ddz6gl2SfJA/vl+wK/Clwy3qpGo6pOqqoDq+oRwMuAL1bVK8Zd1ygluV+SvfrlPYHn\nAmO78mtZhXdVbaX7tOe5wBXAR5bZlQpnAF8BHpPk2iTL7QNRhwHHAM/uL8W6JMnacRc1Qr8AfLGf\n874I+GRVfWHMNS2W5TiFuS9wwcDz96mq+uy4illWlwpK0kqxrEbekrRSGN6S1CDDW5IaZHhLUoMM\nb0lqkOEtSQ0yvLVsJPnj/qtkL+uvET8kyWv7D8TMd9//upB20lLhdd5aFpI8A/hLur+fuiXJg4D7\nABcCT6uqG+e5/3cX0k5aKhx5a7l4KPCj/quA6b9v+TfovlL2S0m+AJDk1CTf6Efok/1tvzek3XOT\nfCXJN5P8Q/9xaJK8pf9jEZcl+fNdfpRSz5G3loU+XL8M3A/4PN1XI5w/88vzk+xdVTf3f7jj88Br\nqurywXZJ9qH7cqy1VXVH/xeZ9gDeCXylqh7X97W6qm7d5Qcr4chby0RV3Ub3F1xeDdwAfCTJsf3m\nwa8oPSrJN4H/C/wicPCQ7g7tb/9KkkuAVwAPA24B/jXJ6Ul+HbhjMY5FWog5/3q81JKq2gacB5zX\n/53PY6c3ASR5BPD7dHPbtyR5L928+DCfq6qjZ96Y5BDgCLopmeP7ZWmXc+StZSHJY5KsGbjpKcA1\nwI+B1f1tq4HbgFuT7Ev3hfrTBttdBByW5FF933smWdNPzTywqj4DvB540mIdjzQfR95aLu4P/E3/\nfdlbgSvpplCOBs5Jcl1VHdFPg2yg+3N5g3/4990z2h0LnJHk3v32P6YL+P+T5D50UzGv2xUHJg3j\nCUtJapDTJpLUIMNbkhpkeEtSgwxvSWqQ4S1JDTK8JalBhrckNcjwlqQG/X8G+YQGaQBB/gAAAABJ\nRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 13 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run the above (or below) cell several times;\n", "you will observe that the distribution differs across sample paths.\n", "Sometimes the state is absorbed into the recurrent class $\\{0\\}$,\n", "while other times it is absorbed into the recurrent class $\\{1, 4\\}$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "init = 2\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(8, 3))\n", "titles = ['Some sample path', 'Another sample path']\n", "titles = [title + ' init={0}'.format(init) for title in titles]\n", "\n", "for ax, title in zip(axes, titles):\n", " draw_histogram(time_series_dist(mc0, init=init), ax=ax, title=title, xlabel='States')\n", "\n", "fig.suptitle('Time series distributions', y=-0.05, fontsize=12)\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAd8AAAD3CAYAAABciF63AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHutJREFUeJzt3XmYXGWZ9/Hvj4SdBKIoSghEEBQUlAEDimAjCpFRUEcG\ng4AsM8I7L+4Loo408qIyzgygOCy+CIoKKI4KDgKitCBLAIewTBIkYiALhC2GJVESc88f5yk4FLV1\nuvqp092/z3XV1afqPPWcu06du++z1TmKCMzMzCyftXodgJmZ2Vjj4mtmZpaZi6+ZmVlmLr5mZmaZ\nufiamZll5uJrZmaWmYvvGCJpvqR9KhDHEZKuH0T7uyXt1e22NjZJWi1p617HMZyc69U36oqvpDdL\nulHSnyQ9Jum3knbtdVwVEemRjaSp6Z/dGi9rEfHaiLhusG0l9Uu6cE2mKelfJf1e0hOS5kg6bE36\nsc5JGpD0uKR1utzn0d3qbwRxrneoV7k+qoqvpInAz4EzgEnAZOAk4C+9jMsAUK8DGKSngHdGxETg\ng8AZkt7Y45hGLUlTgWnAw8ABXex62AqQpPHD1fcI51zvRESMmgewK7C0xXgBXwDmA0uA7wAT07ip\nwGrgCOAB4DHgWOANwJ3AUuAbdf0dBcwGHgeuBLZsMt31gO8Bj6Z+bgFemsYdmfp4AvgD8KHS+/qA\nhcCnKf4pLQbeDewP/D7F+NlS+37gUuDi1N/vgJ1K4/8IvLU0Lz4LzEtxXQJMahJ/LY4TgEdSP4eU\nxv8tcDuwLM27E0vjHkjz9ckU0+4UC/j1wNfSvLsPmN7ie5tfirsf+GH67p4A7gZ2qWu7DzCdYqXr\nmTTt24e4bP0M+ESvl/HR+gC+CFwGfB64vG7cBcA3KVasnwBuBrYujX8TcCvwp5Rbb0yvnwKsAlak\nZeDr6fXVwDEph5YCZ9ZNr2lep/f+E3Av8IcGn8O57lzvbDq9TrqufhiYkBauC9IXMqlu/FEpaaYC\nGwI/Br6bxk1NC85/AOsAb09f6E+ATYHNKQr2Xqn9gamvV1HsQfg8cEOTuI6h+MeyXkqEnYEJadz+\nwCvS8F7A08DOpURYSbHCMA74h/T5vp/i3wFYDmxVWlifAd6b2n8yLezjGiTkR4Eb0+daGzgb+EGL\nhFwJ/GtquxfF2uJ2afxbgNek4R2Bh4AD0/Ot0nxdq9TfESnOo9P8OBZY1OJ7LcfdT/HPdHp675eB\nm5q0PbH2/ZbG/wfFP8VGj1lNpr8+xT/DfXu9jI/WB0Vh+ACwbVo2Xload0Fa7ndNy/X3gIvSuBel\n7+4DKQ/fT/FPflIafy1wVN20VlPk40RgCkWx2y+Na5nX6b1XAZsA6zb4HM5153pny3yvk24YkvjV\nwPnAgrQQ/Yzn1jx/BRxbartdWjDW4rni+/LS+EeBg0rPLwU+koZ/QSmpUx9PA1MaxHQkcAOwYwfx\n/6Q0jb6UcErPJ6QY31BqfxtwQGlhvbE0TmlB2qPBwjq7Npyev7w2LxrEVEvI9UuvXQJ8oclnOB34\n9zRcm6/1CXlv6fkGqc1Lm/RXn5BXl8btACxv0fbCLixT3wGu6PWyPVofwJsp/snWitQs4GOl8ecD\n55aevwOYk4YPA26u6+9G4INp+Frg6Lrxq4E3lZ5fAnwmDbfM6/TevhafxbnuXO/oMaqO+QJExNyI\nODIipgCvpVjbOz2Nfjlwf6n5A8B4YLPSa0tKwysaPN8oDW9FcWxgqaSlFLuFoDjOXO9CirXliyUt\nknRq7XiRpHdIujmdHLaUYu34xaX3PhZpqUjTbxTjRqXnC0vzItLzzRvENBX4SSn+2RS76DZr0BaK\n3fkrSs/vr/UraTdJ10p6WNKfKNb+X9yok5KHSnEuT4MbNWlbr/z5lwPrDeUkj1YkfY0i6f9+OPo3\noNg1eXVEPJme/yi9VtZsmd+cIo/Lnl02k+CFHioNL2dweb2gyecA53ojzvUGRl3xLYuIeyjWZF6b\nXlpMsSDWbEmxEC5h8B6gOGYzqfTYMCJubhDHqoj4UkS8huL41DuBwyWtS7Hr+18o1gQnAVcwtBMW\nptQG0kK6BcXnbhT/9Lr4N4iIB5v0O0nSBqXnWwGL0vAPgJ8CW0TEJhS7tWrLVqN/fLm8YNqSzpb0\nZJPHXXVtTwL2o9gF9VSuoMcSSetT/LN7q6QHJT1IsQv1dZJ26qCLRRTLYll52Rzs8tdJXjft07ne\nMyMu10dV8ZX0KkmfkDQ5PZ8CzABuSk0uAj6eTonfiOIYwsURsXowk0l/zwY+J2mHNK2NJR3UJK4+\nSTtKGkdxQsBK4K8Ux5bXodi9vVrSO4B9BxFLI7tIek9a2/4Y8GeKE1TqnQ18WdKWKcaXSGp3lulJ\nktaWtCfFiRc/Sq9vRLG2/IykacAhPJcMj1DsZtpmSJ9qzTwETJX07D+4iDg2IiY0eexYayfpBIpl\n5+0RsbQHsY8V76ZYAd4eeF16bE9xks7hqU2rAvULYDtJMySNl3QwxaGnn6fxS2i/7Ik1yOuGHTnX\nnesdGlXFl2Jh3w2YKekpiqJ7J8WaNMC3KXYLXUdxcsJy4MOl93ey5hYAEfFT4FSK3UvLgLso1pwa\neRnFwruMYpfPAMXxiSeBj1Cc0fc4xQLws0bT6zDGSO8/OPX3AeC9EfHXBm3PoDgx5GpJT1DMq2kt\n+n6I4kSFxRTz8JiI+H0a90/Al1I//0xxjKgIqNjNdApwg4rfcO5G498gdrrWPJj31v5hPCbptg77\nrzmFYstiXmlt+bOD7MPaOxz4dkQsjIiH02MJcCZwSCpiTb/ziHiMYuvykxSF7VMUPxt5PLU7A3hf\nWvZOp7Fn++8gr9stp85153pHagf3mzeQvk2x5vNweW2hrs3XKU6CWA4cERG3dztQa0/SicArI6Kr\nPxKX1EfxD2RKu7ZWXc7l0cO5PvJ1suV7PsWp3g1J2p9iIdgW+BBwVpdis8EbaT9ut7ycy6OHc32E\na1t8I+J6il0QzRxAcVITETET2ERSs7PobHg12k3Tzb5tBHMujyrO9RGuG5dHm8zzT71fSHHW3Zqc\nQWxDEBEnDVO/AxRnhtvo5lweIZzrI1+3Triq3wXiNSezkcm5bJZBN7Z8F1H6vRnFmvKi+kaSnMRm\nHYqIXhzTcy6bdVmzXO7Glu9lpN/jSdod+FP6qUCjIIb0OPHEE4f9kl8jJY4qxOA4hieOHnIuO46e\nx1CVOIY7l9tu+Uq6iOJi2ptKWkBxAeu1UwKeExFXSNpf0jyKa6AeOYhkN7NMnMtm1dG2+EbEjA7a\nHNedcMxsuDiXzapjRF3hqq+vr9chANWIowoxgOOoV5U4qq4q88lxPJ/jyBdD2ytcdW1CUuSaltlI\nJonozQlXHXEum3WmVS6PqC1fMzOz0cDF18zMLDMXXzMzs8xcfM3MzDJz8TUzM8vMxdfMzCwzF18z\nM7PMXHzNzMwyc/E1MzPLzMXXzMwsMxdfMzOzzFx8zczMMnPxNTMzy8zF18zMLDMXXzMzs8xcfM3M\nzDJz8TUzM8vMxdfMzCwzF18zM7PMXHzNzMwyc/E1MzPLzMXXzMwsMxdfMzOzzNoWX0nTJc2VdK+k\n4xuM31TSlZJmSbpb0hHDEqmZDZnz2awaFBHNR0rjgHuAtwGLgFuBGRExp9SmH1g3Ik6QtGlqv1lE\nrKrrK1pNy8wKkogIDUO/Xcln57JZZ1rlcrst32nAvIiYHxErgYuBA+vaPAhMTMMTgcfqC6+ZVYLz\n2awixrcZPxlYUHq+ENitrs23gF9LWgxMAP6+e+GZWRc5n80qol3x7WTf0ueAWRHRJ2kb4JeSXhcR\nT9Y37O/vf3a4r6+Pvr6+QYRqNjoNDAwwMDCQY1Jdy2fnstkLDSaX2x3z3R3oj4jp6fkJwOqIOLXU\n5grglIi4IT3/FXB8RNxW15ePE5l1YBiP+XYln53LZp0ZyjHf24BtJU2VtA5wMHBZXZu5FCdwIGkz\n4FXAfUML2cyGgfPZrCJa7naOiFWSjgOuAsYB50XEHEnHpPHnAF8Gzpd0B0Ux/0xEPD7McZvZIDmf\nzaqj5W7nrk7Iu6rMOjJcu527xbls1pmh7HY2MzOzLnPxNTMzy8zF18zMLDMXXzMzs8xcfM3MzDJz\n8TUzM8vMxdfMzCwzF18zM7PMXHzNzMwyc/E1MzPLzMXXzMwsMxdfMzOzzFx8zczMMmt5S0GzXpDy\n3tDHd+gxs9xcfK2ichXEyt65z8xGMe92NjMzy8zF18zMLDMXXzMzs8xcfM3MzDJz8TUzM8vMxdfM\nzCwzF18zM7PMXHzNzMwyc/E1MzPLrG3xlTRd0lxJ90o6vkmbPkm3S7pb0kDXozSzrnA+W06Ssj5G\nErW6rq2kccA9wNuARcCtwIyImFNqswlwA7BfRCyUtGlEPNqgr/A1dK0TRRLlu7xk1ZZLSURE1/+T\ndCufncvWKedy81xut+U7DZgXEfMjYiVwMXBgXZtDgB9HxEKARoXXzCrB+WxWEe2K72RgQen5wvRa\n2bbAiyRdK+k2SYd1M0Az6xrns1lFtLurUSfb8GsDfwPsA2wA3CTp5oi4d6jBmVlXOZ/NKqJd8V0E\nTCk9n0Kxtly2AHg0IlYAKyRdB7wOeEGy9vf3Pzvc19dHX1/f4CM2G2UGBgYYGBjIMamu5bNz2eyF\nBpPL7U64Gk9xgsY+wGLgFl54gsargTOB/YB1gZnAwRExu64vn6RhHfFJGsN2wlVX8tm5bJ1yLjfP\n5ZZbvhGxStJxwFXAOOC8iJgj6Zg0/pyImCvpSuBOYDXwrfrCa2a953w2q46WW75dnZDXlq1DXlse\nni3fbnEuW6ecy2v+UyMzMzPrMhdfMzOzzFx8zczMMnPxNTMzy8zF18zMLDMXXzMzs8xcfM3MzDJz\n8TUzM8vMxdfMzCwzF18zM7PMXHzNzMwyc/E1MzPLzMXXzMwsMxdfMzOzzFx8zczMMnPxNTMzy8zF\n18zMLDMXXzMzs8xcfM3MzDJz8TUzM8vMxdfMzCwzF18zM7PMXHzNzMwyc/E1MzPLzMXXzMwss7bF\nV9J0SXMl3Svp+Bbt3iBplaT3djdEM+sW57NZNbQsvpLGAWcC04EdgBmStm/S7lTgSkDDEKeZDZHz\n2aw62m35TgPmRcT8iFgJXAwc2KDdh4FLgUe6HJ+ZdY/z2awi2hXfycCC0vOF6bVnSZpMkcBnpZei\na9GZWTc5n80qYnyb8Z0k3unAZyMiJIkWu6n6+/ufHe7r66Ovr6+D7s1Gt4GBAQYGBnJMqmv57Fw2\ne6HB5LIimuejpN2B/oiYnp6fAKyOiFNLbe7juQTdFFgO/GNEXFbXV7SalllN8T8/17IiqrZcSiIi\nun6stVv57Fy2TjmXm+dyu+I7HrgH2AdYDNwCzIiIOU3anw9cHhH/2WCcE9Y64oQdtuLblXx2Llun\nnMvNc7nlbueIWCXpOOAqYBxwXkTMkXRMGn9O16M1s2HhfDarjpZbvl2dkNeWrUNeWx6eLd9ucS5b\np5zLzXPZV7gyMzPLzMXXzMwss3Y/NTIzGxGKXZz5VG0Xp40sLr5mNorkO75oNhTe7WxmZpaZi6+Z\nmVlmLr5mZmaZufiamZll5uJrZmaWmYuvmZlZZi6+ZmZmmbn4mpmZZebia2ZmltmYvsKVL0dnZma9\nMKaLb8GXozMzs7y829nMzCwzF18zM7PMXHzNzMwyc/E1MzPLzMXXzMwsMxdfMzOzzFx8zczMMnPx\nNTMzy8zF18zMLLOOiq+k6ZLmSrpX0vENxn9A0h2S7pR0g6Sduh+qmQ2Vc9msGtoWX0njgDOB6cAO\nwAxJ29c1uw/YKyJ2Ak4Gzu12oGY2NM5ls+roZMt3GjAvIuZHxErgYuDAcoOIuCkilqWnM4Etuhum\nmXWBc9msIjopvpOBBaXnC9NrzRwNXDGUoMxsWDiXzSqik7sadXzbH0l7A0cBe6xxRGY2XJzLZhXR\nSfFdBEwpPZ9Cscb8POnEjG8B0yNiaaOO+vv7nx3u6+ujr69vEKGajU4DAwMMDAzkmJRz2WwYDSaX\n1e4G75LGA/cA+wCLgVuAGRExp9RmS+DXwKERcXOTfqJqN5OXRM77+Vbt81fVWP9eJBERXb8B9GjO\nZfByU0Vj/Ttplcttt3wjYpWk44CrgHHAeRExR9Ixafw5wBeBScBZxcxmZURM69YHMLOhcy6bVUfb\nLd+uTaiCa8tjfa2sqsb69zJcW77dUsVcBi83VTTWv5NWuewrXJmZmWXm4mtmZpaZi6+ZmVlmLr5m\nZmaZufiamZll5uJrZmaWmYuvmZlZZi6+ZmZmmbn4mpmZZebia2ZmlpmLr5mZWWYuvmZmZpm5+JqZ\nmWXm4mtmZpaZi6+ZmVlmLr5mZmaZufiamZll5uJrZmaWmYuvmZlZZi6+ZmZmmbn4mpmZZebia2Zm\nlpmLr5mZWWYuvmZmZpm5+JqZmWXWtvhKmi5prqR7JR3fpM3X0/g7JO3c/TALAwMDw9X1IA30OoDK\nzIuqxFGF7wSqND8aq0o+V2c+DfQ6AKA686MqcVThexnuedGy+EoaB5wJTAd2AGZI2r6uzf7AKyNi\nW+BDwFnDFOuoXTAkDfqx9957r9H7JHU19tH6nayp6syPF6pSPldnPg30OgCgOvOjKnFU4XvpafEF\npgHzImJ+RKwELgYOrGtzAPAdgIiYCWwiabOuRzrqxSAfJ67BeyLbp7FKcj4PszVdIT7ppJMqsTJt\n+bQrvpOBBaXnC9Nr7dpsMfTQzKzLnM9ZrMlKsVemx5yIaPoA/g74Vun5ocA36tpcDuxRen4N8DcN\n+lqTJcsPP8bko1VerumDLuVzr+eNH36MpEezfBxPa4uAKaXnUyjWhFu12SK99jwR4f0jZr3VlXx2\nLpsNXbvdzrcB20qaKmkd4GDgsro2lwGHA0jaHfhTRCzpeqRmNlTOZ7OKaLnlGxGrJB0HXAWMA86L\niDmSjknjz4mIKyTtL2ke8DRw5LBHbWaD5nw2qw6lYzhmZmaWyYi4wlUnFwbIFMe3JS2RdFcPY5gi\n6VpJ/yPpbkkf6VEc60maKWmWpNmSvtKLOFIs4yTdLunyXsWQ4pgv6c4Uyy29jKXKqpDPVcjlFEfP\n87lKuZzi6Xk+58jlym/5pgsD3AO8jeLEj1uBGRExpwex7Ak8BXw3InbMPf0Uw8uAl0XELEkbAb8D\n3t2j+bFBRCyXNB74LfCpiPhtD+L4BLALMCEiDsg9/VIcfwR2iYjHexVD1VUln6uQyymOSuRzVXI5\nxdLzfM6RyyNhy7eTCwNkERHXA0t7Me1SDA9FxKw0/BQwB9i8R7EsT4PrUBxDzF50JG0B7A/8f6AK\nZ+FWIYYqq0Q+VyGXUxyVyOcq5DJULp+Hdfojofh2cmGAMUnSVGBnYGaPpr+WpFnAEuDaiJjdgzBO\nAz4NrO7BtOsFcI2k2yT9Y6+DqSjncxO9zOeK5DJUJ5+HPZdHQvGt9n7xHkm7qC4FPprWmLOLiNUR\n8XqK34LuJakv5/QlvRN4OCJup/dryVBcnGJn4B3A/027Nu35nM8N9Dqfe53LULl8HvZcHgnFt5ML\nA4wpktYGfgx8LyJ+2ut4ImIZ8F/Arpkn/SbggHR85iLgrZK+mzmGZ0XEg+nvI8BPKHax2vM5n+tU\nKZ97mMtQoXzOkcsjofh2cmGAMUPFldTPA2ZHxOk9jGNTSZuk4fWBtwO354whIj4XEVMi4hXA+4Ff\nR8ThOWOokbSBpAlpeENgX6CnZ9JWlPO5pAr5XIVchurkc65crnzxjYhVQO3CALOBS3pxZi+ApIuA\nG4HtJC2Q1IsLEOxBcU3evdNp8LdLmt6DOF4O/DodJ5oJXB4Rv+pBHGW93KW5GXB9aX78PCKu7mE8\nlVSVfK5ILkM18rmKuQy9y+csuVz5nxqZmZmNNpXf8jUzMxttXHzNzMwyc/E1MzPLzMXXzMwsMxdf\nMzOzzFx8zczMMnPxHQUkfT7djuyO9DvBaZI+mn4w3+69H+uknZnl4XweG/w73xFO0huBfwPeEhEr\nJb0IWA+4Adg1Ih5r8/4/dtLOzIaf83ns8JbvyPcy4NF0ezbS/SffR3Fbsmsl/QpA0lmSbk1r1P3p\ntY80aLevpBsl/U7SD9Pl1ZD0VRU3/L5D0teyf0qzscH5PEZ4y3eES8n0W2AD4BqKy/VdV38zaEmT\nImKpipuZXwN8OCLuLreTtCnFBd6nR8QKScdT3N/zm8CNEfHq1NfEiHgi+4c1G+Wcz2OHt3xHuIh4\nGtgF+BDwCHCJpCPS6PJtuQ6W9Dvgv4HXADs06G739PqNkm4HDge2BJYBf5Z0nqT3ACuG47OYjXXO\n57FjfK8DsKGLiNXAb4DfSLoLOKI2CkDSK4BPUhwLWibpfIrjSI38MiIOqX9R0jRgH4pdYMelYTPr\nMufz2OAt3xFO0naSti29tDMwH3gSmJhemwg8DTwhaTOKG0TXlNvNBPaQtE3qe0NJ26ZdYZtExC+A\nTwCvG67PYzaWOZ/HDm/5jnwbAd9QcT/OVcC9FLusDgGulLQoIvZJu53mAgsojinVnFvX7gjgIknr\npvGfp0jon0laj2LX18dzfDCzMcj5PEb4hCszM7PMvNvZzMwsMxdfMzOzzFx8zczMMnPxNTMzy8zF\n18zMLDMXXzMzs8xcfM3MzDJz8TUzM8vMxdfMzCwzF18b0dL9TPfqdRyDIekKSYcN8zT6JV2YhreU\n9KQktXtfh32fJekLabhP0oJu9Jv621PS3G71Z1ZVvrazVZqkp0h3cwE2BP4M/DU9/1BEvLYngQ1B\nROyfYzKl6T0ATGj3hnQd4KMjYs+WHUf8nyFH99w0VwOvjIj7Ut/XA6/uVv9mVeXia5UWERvVhtON\nwo+OiF/3MKQ1VtvyjBF8QXVJa6Vb3nW12y73Z1Z53u1sI5qk+ZLemob7Jf1I0oWSnpB0Z7qF2gmS\nlki6X9LbS+/dON1QfLGkhZJOltQwJyRNk3SbpGWSHpL0b6Vxu0u6UdJSSbMkvaU0bkDS/5N0A/AU\nsHV67ehSm6MkzZb0uKQrJW1ZGndain1Z+jyvaRLfKyT9Jn3uq4FNS+OmSlpd+2ySjpD0h9T2PkmH\nSHo1cDbwxrSL+vHU9oK0m/mKtBdi7/TayXXTP0HSI5L+KOmQ0uv1n/UISden4evSy3ekaR5Uvxtb\n0vapj6XpEMO7SuMukPRNST9Pn+VmSVsPdt6Z9YKLr4109VuR7wS+C0wCbgd+mV7fHDgZOKfU9gLg\nGWAbivum7gv8Q5PpnAGcFhEbA1sDPwSQNBn4OfCliJgEfAr4saQXl957aOp3AnB/irl2Y/QDgROA\n91AUzOuBi9K4/YA9gW3TdA8CHmsS3w+AW4EXp8/5wQbzBhX3cj0DmB4RE4E3ArMiYi5wDHBTREyI\niBeV3jYDODnthfhtOf7kZWm6m6fpnqvn7klb3/ZZEVE7Vr9TmuaP6mJdG7gcuBJ4CfBh4PuStis1\nOxjop/i+5wGnpPcOZt6ZZefia6PNdRHxy4j4K3ApRVH4anp+CTBV0kQ9dxPyj0fEioh4BDgdeH+T\nfp8BtpW0aUQsj4iZ6fVDgSsi4kqAiLgGuA342zQ+gAsiYk5ErI6IVXX9Hgt8JSLuSbtzvwK8Pm39\nPkNRsLdPu3vviYiH6gNLbXcF/jkiVqbjppfTfHfuamBHSetHxJKImF3rqkHbAH4aETelz/eXJm1r\n074O+C+KojhUuwMbRsRXI2JVRFxLsaIzo9TmPyPitvT9fh94fXp9JR3MO7NecfG10ebh0vAK4NHS\nMdYV6e9GwFbA2sCDaZfmUordri9p0u/RwHbAHEm3SKoV162Ag2p9pH72oNgarGl1NvBWwBml99a2\nzjZPxeZM4JvAEknnSGp04tTmwNKIWFF67f5GE4uIpykK47HA4rTL9lUt4msXP02m/fI27+nE5g2m\nfX96HYoVgyWlcSsovlvSeQGdzDuznnDxtbFqAfAX4MURMSk9No6IHRs1joh5EXFIRLwEOBW4VNIG\nwAPAhaU+JqVdqP9SfnuLOB6gOGu7/P4NI+LmNN1vRMSuwA4Uxf/TDfp4EJiU4qnZqtl0I+LqiNiX\nYgVhLvCtDuJ8QTel4UbTXpyGn6Y4S72mvFLSzmJgivS8n0htBSzqKMDO5p1ZT7j42pgUEQ8CVwP/\nLmmCpLUkbaMmvxmWdKik2lbxMori81fge8C7JO0raZyk9dJJQ5PLb28RytnA5yTtkKazsaSD0vCu\nknZLxz6X8/yfWZU/y/0Uu7pPkrS2pDdTHPtu9DleKunAdOx3JUVxrPW5BNgiTa9V7Grwem3ae1Ls\ncq8dv50FvFfS+pJeSbEHoWwJxTH3RmZSfO7PpL770ue6uEVstc/Z0bwz6xUXXxtNGp3c0+r54cA6\nwGzgcYqC0WzLbD/gbklPAqcB74+Iv0TEQuBA4HMUu7wfAD7J8wtD0y3KiPgpxZb0xZKWAXelaQFM\nBM5Nsc0HHgW+1qSrQ4DdUtsvAt9p8rnXAj5OsfX4GMVJSbXf7f4K+B/gIUkPl97XaB6WX3sQWEqx\npXohcExE/D6NO43i2PUS4HyKlZXye/uB76Td7u8r9x0RzwDvojg2/wjFbuTDSn23+r4HM+/MstMI\n/smhmZnZiOQtXzMzs8xcfM3MzDJz8TUzM8vMxdfMzCwzF18zM7PMXHzNzMwyc/E1MzPLzMXXzMws\nMxdfMzOzzP4Xoxlo4F5BXK0AAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 14 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us repeat the simulation for many times (say, 100 times)\n", "and obtain the distribution of visits to each states\n", "at a given time period `T`.\n", "That is, we want to simulate the marginal distribution at time `T`." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def cross_sectional_dist(mc, init, T, num_reps=100):\n", " \"\"\"\n", " Return a distribution of visits at time T across num_reps sample paths\n", " of mc with an initial state init.\n", " \n", " \"\"\"\n", " x = np.empty(num_reps, dtype=int)\n", " for i in range(num_reps):\n", " x[i] = mc.simulate(init=init, sample_size=T+1)[-1]\n", " bins = np.arange(mc.n+1)\n", " hist, bin_edges = np.histogram(x, bins=bins)\n", " dist = hist/len(x)\n", " return dist" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 15 }, { "cell_type": "code", "collapsed": false, "input": [ "init = 2\n", "T = 100\n", "draw_histogram(cross_sectional_dist(mc0, init=init, T=T),\n", " title='Empirical marginal distribution ' + \\\n", " 'at T={0} with init={1}'.format(T, init))\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAAEKCAYAAADdBdT9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGGpJREFUeJzt3HmUXGWdxvHvQ0JYZBcFCWFwARVxQRRRBFpwCYwDqMfR\nCCLICOPIOHNwQRwHWlBR5xwXBg8yGlFRiYwKBkWiAsW+xUNASAJEyZiFfQ0EJTG/+eN9O16qa+vu\nqlS/nedzTp1U1b393t/dnnrrvbeiiMDMzMqyQb8LMDOzkXN4m5kVyOFtZlYgh7eZWYEc3mZmBXJ4\nm5kVaL0Pb0krJO3cYvpZkj4zxmUMSFoyljZ6RdLhkuZ0qa3Fkg7scN6jJF1Ved1yP4ywjpMkfSs/\n31nSGkldOdYl7ZRrVTfas/Yk3SZpvxbTa5KO6bCtfSUt7Pa8/TAuwzuHwMp8kgw9zujFsiJi84hY\n3GL6hyPic71Y9ngQET+MiLd1q7n8GE0dLfcDdP4hGBGnR8SHRlNHg2UulnRApe0/5VrX6Q8kJA1K\nOrfF9Ccq58qauvNnxgiWs7ukOZIekLSmwfRtJF2Ql7e4vm1JB0paKOlJSZdJ2mlkazpcROweEVfm\n9htth46Pu4i4KiJeMpp564+FTknaVdLPJd0v6SFJl0jadaTt1BuX4U3aEW/PJ8nQ46Pruohu9db6\nRVm/61iXJE3qcpMBjPttGBGbDZ0rwP/xzPPnvBE09TQwC2jWk/0G8GfgucDhwFmSdgOQtC3wU+A/\ngK2BucCPR7VC49Noj4UtgQuBXYHtgBuBn4+9mohx9wDuBg5oMu0o4BrgK8AjwCLgDcDRwJ+A+4Aj\nK/N/F/gm8GvgcaAG7FSZvgZ4QWXes4CLgSeAA/N7p1XmPxSYBzyWl/22/P7RwPy8jD8Ax1b+ZgBY\n0mJ91wAfBu7Kf38q8ELgOuBR0sm0YZ53K+AXwP3Aw8BFwNRKWzXgc3kbrQReALwVuCO39Q3gCuCY\nyva8qq6W44A78/Y9szLthcBlwIPAA8APgC073G/PBmbn7XYDcFqD5Q7th4OB2/O2WAqcAGwKPAX8\nFViRpz0PGAR+Apyb2z4mv3dubmvn3PaHgGXAcuBjdcfHaY32VW7zr3k7rgA+XmlvgzzPDnm9Hsr7\n758qbQ0C5wPfy/XeBuzZ4jj4OukYfowUfG/M708H/kIK1hXAzaM9f0ZwDr4IWFP33rNyHS+qvPc9\n4PT8/Fjg6sq0TfO227VB+28Cbq28/g1wY+X1VcAh+fli0rnYcDsAl5POmavzdp4DPLvJeq3dv5W2\nPwbcwt/OtY06ORbGsG23ycfQ1mPaR2P541498sF3YJNpRwGrgA+QPgVPI53g/w1sCLwl78BN8/zf\nza/fCEwBvkbz0Phu3oGvz683As4BTs2v98rTD8yvdwBenJ8fDDw/P98PeBLYo9EB02Cd1gAXAJsB\nu+UD9DJSUGxBCrIjKzv+HcDGef7zgQsqbdXyAflS0jer55DC4LD8+qP54P9gZXvWb4/ZebnTSB8S\nQx9QLySdRBsC25I+BL5at9+ahfes/NgEeFneZ1c22Q/3APvk51tWtuP+9duRFJBP87cTfWPgFIaH\n9w/zsnfP6zS0D9fu3yYn9zPWieHhfSVwJunYemVu+02V2p4ihY6ALwDXtTgODif1WDcgfWDdA0zJ\n004Bvj+C8+eAuvfeR/owbvR4GNixbv5G4b0H8GTdeycAs/PzrwPfqJt+K/DOBjVukrfNNvl4ug9Y\nQvqA2IQUklvXr0+j7UA65hflmjcmhfnpTbZNo/17PbB93vbzgeM6ORbye4+22K6fbFLDYcCysWRk\nRIzbYRMBF0p6pPKofo27OyK+F2lLnE8K0VMjYlVE/IZ0Mr+oMv8vIuLqiHia9JXu9ZKmNln2hRFx\nHUBE/KVu2jHAzIi4NE9fHhF35OcXR8Td+fmVpJ7+viNY5y9HxBMRMR/4PfCriFgcEY8DvyKdOETE\nwxFxQUT8OSKeIAXC/pV2AvhuRCyIiDXAQcBtEXFhRKyJiDOAe9vU8sWIeDwilpBOhFflZf8hIi7N\n2/lB4Kt1y24oD2W8Ezg5Ip6KiNtJPbZmX0GfBl4maYuIeCwibh5qqsn810bE7Fzjn5vM99m87NtI\ngV0dqx3VsIikaaRvfSdGxNMRcQvwbeDIymxXRcQl+Vj9ASngG4p0/eGRvJ++Quo8vLhS46iHbyLi\nRxGxdZPHNhGxtINmNiN1hKpWAJu3mP54fr++nqeAm0jHz56kb7PXkDpZewN3RcQjDWpotB0C+E5E\nLMr7/3zyMduhMyLi3ry8i0bytxGxVYvt+uVhxUs7kj7sTxhBfQ2N1/AO4NC6DTGzMv2+yvOnACLi\ngbr3hg6YIPXyyPM9Sepp7NBkua0uiO1IGhIZRtJBkq7PFyQeIfXEn92irXr161T/erO8nE0lnZ0v\nnjxG6v1uWTe2XV2HHaisf9buRK2G+8rKsreTNEvS0rzsc+lsHZ8DTK6r608t5n8XafstzncS7N2m\n/U6Cp37Zjfb/SO0APJyPqWrb1Y5BdT+uBDZudi1F0sclzZf0aD6GtiR9wxkvniB9I6vakr8F9oom\n01c0ae8KUu923/z8ClKY70fqTY9E9Zitnv+9/tuOSXoOqVP3jYgY87WA8Rre3STS1//0QtqM9FVt\n+SjaWsIze/RDbW5EulDzZeC5EbE1ady8Fxe6Pka68LFXRGxJOtjreyNReb6c9KEzVKuqrzs01N4X\nSON+u+dlv5/OjqEHgNVA9c6DpnchRMTciDiMFPoXknpS1Trqa6t/v9F89ctelp8/SRqbHbJ9B20N\nWQ5sk4+patudfJg8g6R9gU8A7x7qzZGGu4b2a6s6Omn/8Lq7t6qPx3OPsJ07gcmSqufAK0nDeuR/\n136zkPQs0lDb7TR2BWnseyish8J8//y8kTFthzEatuy6u3zqH5+qzLc1KbgvjIjTu1HMeA7vbgbf\nwZL2kTSFNEZ+XUQsazBfo2VWg3EmcLSkAyRtIGmqpBeTxjunkC7krZF0EOki4VioyfPNSL2DxyRt\nQxoDbPW3vwReLulQSZOBjzA8oDqtYzNS2D2eh50+0UkDEfFX4GfAoKRN8t0JH2i4MGnDHDRb5r9b\nQfrAgNSLfbakau+u2T6r95m87JeRxvmHej7zSMfH1pK2B/697u/uIwVQo/VaAlwLnC5pI0mvAD5I\nGh4Zqc1JH3APSpoi6WSe2Yu9F9h5tHcP5SGZzZs8tqgOm0jamHQ8k9dro9zGk6T9eGr+BvhG4B9I\n38AgXbfZXdI7cxunAPMi4s4mZV1LGhZ6Leli5Xzg74DXka4lNNJsO6yLO4KGHQtRucunweOLAPl4\nnUO6mPvpbhUznsP7orpPsZ/m9zvtaVWn/Yh0ID1EGjs+osnfNms7ACLiJtJdJV8lXaioke5cWUG6\nEHg+aUhmBsNvBWpXY6v3qnV9jXRB50HSwf+rJjWTa34IeDfpW8GDpAuZc0kXRevbblRL9fVngVeT\neoQXkb5tdNoTOp4U/vcC38mPZss9Arg7D80cS7qQR0QsBM4D/ijpYUnPa1B/s3W6gnRR67fAf0XE\nb/O0c0l3GiwGLiFdVK3+7emk4H9E0gmV9obMIF3EXE4KtpMj4rImddT/bdUl+XFnruUpnjm09L/5\n34ckzW3Sxpgp/VBqJenOmMh1LKjM8i+k4+9+0ofUP0fEAoB8HeRdwOdJ58FrgPc2W1ZErAR+B9we\nEavz29cCi3NbjTTbDu3OY+qmt5rW7LhsdCx04h2kbXH0KL7tNKV0HaXFDNJ3gL8H7o+IlzeZ5wzS\nhbGVwFGVC0x9J+kcYGlE/Ge/axkP8njrEuB9EdHsq6mZjXOd9LzPId3q1JCkg0n3fe5C6iWd1aXa\numXc/8Ci1yS9VdJW+evv0Ne26/tZk5mNTdvwjoirSPcsNnMI6bYvIuIGYCtJ23WnvK5o9xVqffB6\n0pDBA6RvUYc1uA3SzAoyuQttTOWZt2EtJd3NcF/j2detiDi63zX0W0R8ljRebWYTRLcuWDa6ad7M\nzHqkGz3vZVTuoyb1uofdhifJgW5mNgoRMezaXTd63rPJPwfOv4R7NCIaDpnEGH/LP5LHKaecsk6X\nt64fE3n9JvK6ef3Kf6zr9Wumbc9b0nmkXzxtq/R/KZ9C+o9kiIizI+JiSQdLWkT6Acd6P8ZsZtZr\nbcM7Itr+R+4RcXx3yjEzs06M519YjsnAwEC/S+ipibx+E3ndwOtXuvGyfm1/Ydm1BUmxrpZlZjZR\nSCJ6dMHSzMzWMYe3mVmBHN5mZgVyeJuZFcjhbWZWIIe3mVmBHN5mZgVyeJuZFcjhbWZWIIe3mVmB\nHN5mZgVyeJuZFcjhbWZWIIe3mVmBHN5mZgVyeJuZFcjhbWZWIIe3mVmBHN5mZgVyeJuZFcjhbWZW\nIIe3mVmBHN5mZgVyeJuZFcjhbWZWIIe3mVmBHN5mZgVyeJuZFcjhbWZWIIe3mVmBHN5mZgVyeJuZ\nFcjhbWZWIIe3mVmB2oa3pOmSFkq6S9KJDaZvK+kSSfMk3SbpqJ5UamZmaykimk+UJgF3AG8GlgE3\nATMiYkFlnkFgo4g4SdK2ef7tImJ1XVvRallmZjacJCJC9e+363nvBSyKiMURsQqYBRxaN889wBb5\n+RbAQ/XBbWZm3TW5zfSpwJLK66XA6+rm+RZwmaTlwObAP3avPDMza6RdeHcyzvFpYF5EDEh6IfAb\nSa+MiBX1Mw4ODq59PjAwwMDAwAhKNTOb+Gq1GrVare187ca89wYGI2J6fn0SsCYivlSZ52Lg8xFx\nTX59KXBiRMyta8tj3mZmIzTaMe+5wC6SdpY0BXgPMLtunoWkC5pI2g54MfDHsZdsZmbNtBw2iYjV\nko4H5gCTgJkRsUDScXn62cAXgHMk3UL6MPhkRDzc47rNzNZrLYdNurogD5uYmY3YaIdNzMxsHHJ4\nm5kVyOFtZlYgh7eZWYEc3mZmBXJ4m5kVyOFtZlYgh7eZWYEc3mZmBXJ4m5kVyOFtZlYgh7eZWYEc\n3mZmBXJ4m5kVyOFtZlYgh7eZWYEc3mZmBXJ4m5kVyOFtZlYgh7eZWYEc3mZmBXJ4m5kVyOFtZlYg\nh7eZWYEc3mZmBXJ4m5kVyOFtZlYgh7eZWYEc3mZmBXJ4m5kVyOFtZlYgh7eZWYEc3mZmBXJ4m5kV\nyOFtZlagtuEtabqkhZLuknRik3kGJN0s6TZJta5XaWZmz6CIaD5RmgTcAbwZWAbcBMyIiAWVebYC\nrgHeFhFLJW0bEQ82aCtaLcvMzIaTRESo/v12Pe+9gEURsTgiVgGzgEPr5nkf8NOIWArQKLjNzKy7\n2oX3VGBJ5fXS/F7VLsA2ki6XNFfS+7tZoJmZDTe5zfROxjk2BF4NHAhsClwn6fqIuGusxZmZWWPt\nwnsZMK3yehqp9121BHgwIp4CnpJ0JfBKYFh4Dw4Orn0+MDDAwMDAyCs2M5vAarUatVqt7XztLlhO\nJl2wPBBYDtzI8AuWLwHOBN4GbATcALwnIubXteULlmZmI9TsgmXLnndErJZ0PDAHmATMjIgFko7L\n08+OiIWSLgFuBdYA36oPbjMz666WPe+uLsg9bzOzERvtrYJmZjYOObzNzArk8DYzK5DD28ysQA5v\nM7MCObzNzArk8DYzK5DD28ysQA5vM7MCObzNzArk8DYzK5DD28ysQA5vM7MCObzNzArk8DYzK5DD\n28ysQA5vM7MCObzNzArk8DYzK5DD28ysQA5vM7MCObzNzArk8DYzK5DD28ysQA5vM7MCObzNzArk\n8DYzK5DD28ysQA5vM7MCObzNzArk8DYzK5DD28ysQA5vM7MCObzNzAo0ud8FSOp3CSMWEf0uwczW\nc30P76SkMCzvw8bMJp62wyaSpktaKOkuSSe2mO+1klZLemd3SzQzs3otw1vSJOBMYDqwGzBD0kub\nzPcl4BLcNTUz67l2Pe+9gEURsTgiVgGzgEMbzPevwE+AB7pcn5mZNdAuvKcCSyqvl+b31pI0lRTo\nZ+W3ShrANjMrUrvw7iSIvwZ8KtItGMLDJmZmPdfubpNlwLTK62mk3nfVnsCsfMvftsBBklZFxOz6\nxgYHB9c+HxgYYGBgYOQVm5lNYLVajVqt1nY+tbpnWdJk4A7gQGA5cCMwIyIWNJn/HOCiiPhZg2nR\naFkp9EsaaZHv8zazdUYSETFsRKNlzzsiVks6HpgDTAJmRsQCScfl6Wf3pFozM2upZc+7qwtyz9vM\nbMSa9bz9f5uYmRXI4W1mVqBx8n+bmNm65v8UrmwOb7P1WklhWN6HTS952MTMrEAObzOzAjm8zcwK\n5PA2MyuQw9vMrEAObzOzAjm8zcwK5PA2MyuQw9vMrEAObzOzAjm8zcwK5PA2MyuQw9vMrEAObzOz\nAjm8zcwK5PA2MyuQw9vMrEAObzOzAjm8zcwK5PA2MyuQw9vMrEAObzOzAjm8zcwK5PA2MyuQw9vM\nrEAObzOzAjm8zcwK5PA2MyuQw9vMrEAObzOzAjm8zcwK5PA2MytQR+EtabqkhZLuknRig+mHS7pF\n0q2SrpH0iu6XamZmQ9qGt6RJwJnAdGA3YIakl9bN9kdgv4h4BXAa8D/dLtTMzP6mk573XsCiiFgc\nEauAWcCh1Rki4rqIeCy/vAHYsbtlmplZVSfhPRVYUnm9NL/XzDHAxWMpyszMWpvcwTzRaWOS3gR8\nENin0fTBwcG1zwcGBhgYGOi06SJJ6ncJoxLR8S43sy6r1WrUarW286ndiSppb2AwIqbn1ycBayLi\nS3XzvQL4GTA9IhY1aCcaLSsFXElhoY7Drbx1g5Gsn5WtvONz/Tw2JRERw3qCnQybzAV2kbSzpCnA\ne4DZdY3vRAruIxoFt5mZdVfbYZOIWC3peGAOMAmYGRELJB2Xp58NnAxsDZyVhwpWRcRevSvbzGz9\n1nbYpGsL8rBJIdbPr6bro/KOz/Xz2BzLsImZmY0zDm8zswI5vM3MCuTwNjMrkMPbzKxADm8zswI5\nvM3MCuTwNjMrkMPbzKxADm8zswI5vM3MCuTwNjMrkMPbzKxADm8zswI5vM3MCuTwNjMrkMPbzKxA\nDm8zswI5vM3MCuTwNjMrkMPbzKxADm8zswI5vM3MCuTwNjMrkMPbzKxADm8zswI5vM3MCuTwNjMr\nkMPbzKxADm8zswI5vM3MCuTwNjMrkMPbzKxADm8zswI5vM3MCtQ2vCVNl7RQ0l2STmwyzxl5+i2S\n9uh+mWZmVtUyvCVNAs4EpgO7ATMkvbRunoOBF0XELsCxwFk9qnWEav0uoMdq/S6gZ2q1Wr9L6KmJ\nvn4T+diE8bP/2vW89wIWRcTiiFgFzAIOrZvnEOB7ABFxA7CVpO26XumI1fpdQI/V+l1Az4yXk6NX\nJvr6TeRjE8bP/msX3lOBJZXXS/N77ebZceylmZlZM+3COzpsR6P8OzMzGwVFNM9ZSXsDgxExPb8+\nCVgTEV+qzPNNoBYRs/LrhcD+EXFfXVsOdDOzUYiI+g4yk9v8zVxgF0k7A8uB9wAz6uaZDRwPzMph\n/2h9cDdbuJmZjU7L8I6I1ZKOB+YAk4CZEbFA0nF5+tkRcbGkgyUtAp4Eju551WZm67mWwyZmZjY+\nTbhfWHbyo6JSSfqOpPsk/b7ftfSCpGmSLpd0u6TbJH203zV1k6SNJd0gaZ6k+ZJO73dN3SZpkqSb\nJV3U71p6QdJiSbfmdbyxr7VMpJ53/lHRHcCbgWXATcCMiFjQ18K6RNK+wBPA9yPi5f2up9skbQ9s\nHxHzJG0G/A44bKLsPwBJm0bESkmTgauBj0fE1f2uq1sknQDsCWweEYf0u55uk3Q3sGdEPNzvWiZa\nz7uTHxUVKyKuAh7pdx29EhH3RsS8/PwJYAGwQ3+r6q6IWJmfTiFdR+p7CHSLpB2Bg4FvM/z24Ylk\nXKzbRAvvTn5UZAXIdzjtAdzQ30q6S9IGkuYB9wGXR8T8ftfURV8FPgGs6XchPRTAbyXNlfShfhYy\n0cJ74owBrcfykMlPgH/LPfAJIyLWRMSrSL9C3k/SQJ9L6gpJbwfuj4ibGSc90x7ZJyL2AA4CPpKH\nMvtiooX3MmBa5fU0Uu/bCiFpQ+CnwA8i4sJ+19MrEfEY8EvgNf2upUveABySx4TPAw6Q9P0+19R1\nEXFP/vcB4ALSUG1fTLTwXvujIklTSD8qmt3nmqxDkgTMBOZHxNf6XU+3SdpW0lb5+SbAW4Cb+1tV\nd0TEpyNiWkQ8H3gvcFlEHNnvurpJ0qaSNs/PnwW8FejbnV8TKrwjYjXp155zgPnAjyfYnQrnAdcC\nu0paImmi/SBqH+AI4E35VqybJU3vd1Fd9DzgsjzmfQNwUURc2ueaemUiDmFuB1xV2X+/iIhf96uY\nCXWroJnZ+mJC9bzNzNYXDm8zswI5vM3MCuTwNjMrkMPbzKxADm8zswI5vM3MCuTwNjMr0P8Db9Pi\na16VJOQAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 46 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Observe that the distribution is close to a convex combination of\n", "the stationary distributions `[1, 0, 0, 0, 0, 0]` and `[0, 1/3, 0, 0, 2/3, 0]`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us simulate with the remaining states, 3, 4, and 5." ] }, { "cell_type": "code", "collapsed": false, "input": [ "inits = [3, 4, 5]\n", "\n", "fig, axes = plt.subplots(1, 3, figsize=(12, 3))\n", "\n", "for init, ax in zip(inits, axes):\n", " draw_histogram(time_series_dist(mc0, init=init), ax=ax,\n", " title='Initial state = {0}'.format(init),\n", " xlabel='States')\n", "\n", "fig.suptitle('Time series distributions', y=-0.05, fontsize=12)\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAr4AAAD3CAYAAAD7T1f4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHrxJREFUeJzt3XucZGV95/HP1wHjhYsoBuVq1DFq4i2yIy4xtpLggBdy\ncxUWXRISSV5Lkk3UsBoTx5BEzWXRaJagS7xgVpJokgVDwAu0IChKAioyw0IUYbgMIIhyiQ7OL3/U\naT3TdnVXz1TXpc/n/XrVa6rqPHWe51TPt+pXTz19OlWFJEmStNo9YNwDkCRJkkbBwleSJEmdYOEr\nSZKkTrDwlSRJUidY+EqSJKkTLHwlSZLUCRa+I5LknCSvWGT7qUneMOC+ZpMcP7zRSWozr9J0MbMa\nlIXvTkhyXZLDBmlbVUdW1RnN445LctG87b9aVX8wYNfVXJYz1sck2ZZk4J95c3zPX04/w5DkgiS3\nJvlGko1JfnnUY9DqY15XVpLnNmM+eVxj0OpiZldG0++9Sb7ZXM4d9RjGaZdxD2DKLTscEyDLaFvL\nbD8svw5sqqqtSdYBFya5sKquHsNYtHqY1xWSZFfg7cBnmL7nWJPLzK6MAl5UVeePoe+xc8Z3SJpP\nmJ9K8idJ7kjy5STrW9tnkxyf5InAXwLPbj5p3dFsf+/cTEmSvZJ8pJn1vCPJ2Un2G3Ac65JcluSu\nJLck+dNm04XNv19v+n1WksclOT/J7UluS/KBJHs2+zkDOBA4u2n/mub+Q5JckuTOJFckee4wnr+2\nqvpiVW1t3XU38I1h96PuMq9D92rgXOBqxlR8a3Uzs0PX2Zxa+A7XOmAT8Ajgj4HTW9sKqKraBJwA\nfLqqdq+qh7e3N9fTPPbA5nIf8M4Bx/B24JSq2hN4LPB3zf3Paf7ds+n30ub2HwKPBp4EHABsoDfQ\nVwDX0/tUuHtV/WnzwvAR4Perai/gNcCHk+y90ECaF5Y7+1zOWuwgmsfeB8wCv1hVNw94/NKgzGvL\njuY1yUHALwAn0+E3U42EmW3ZmfdY4K+bwv+8JE8d8NhXBQvf4fpqVZ1eVQW8H3h0kh9coF2/N4cA\nVNUdVfUPVfXvVXU38EfAoJ/6vg2sTbJ3Vd3bCt/39VlV/1ZVn6iqrVV1O3DKEv0cC5xTVec2j/84\ncBlw5EKNq+pFVbVXn8tLFjuIqnoRsBvwSuC9SQ5c4ril5TKv2+9/R/P658AbquoepvOraU0PM7v9\n/nc0s8cABzWXC4Dz5maiu8DCd7humbtSVfc2V3db7k6SPCTJaektQL8L+CSwZ5JBZlOOB54AbEzy\n2SQvXKSffZKcmWRz088Z9D5J93MQ8NL2p0rgUOBRAx/cMlTVd6rqQ8ClwM+sRB/qNPO6k5K8GNit\nquZmvYKzvlo5ZnYIqurTVfWtqrqvqt4CfJ3vzVivev5y23j0mxGZu//V9IK1rqpuTfJ04F/pvaEs\nOptSVdfS+zRHkp8DPpTk4X0e90fAd4AfraqvJ/lp4B2LjPN64IyqetViY5iT5J+BH++z+cKq6vuC\nMc+uwD0DtpWGzbz2z+vzgYOTzC1F2hP4TpIfrSo/rGpczOzy3mM79S2NM77jsQXYP73fhJ7TninZ\njd6ao7uaQL1xgX0s+Mk0ybFJHtncvIvef+htwG3Nv49rNd+NXkH5jWZt0WsXGGe7/QeAFyc5PMma\nJA9KMpM+vxRQVUc0a5cWuiwYyCQ/nOSIJA9OsmuSY4GDgY8u1F4aAfPa/w30d4G1wNOApwNnAe+i\nt+ZXGhcz2/899oAkhyZ5YLP/19Kbhb54ofarkYXv8Cy0tq3fp6hPAF8Cbkly6wKPfxvwYOB24BLg\nn5ex7xcAVyb5Jr31RC9vvtK4l94i+4vT+y3WdcCbgB+jF96zgQ/P2++bgTc0X7n8VlVtBo4CXg/c\nSu/T6asZ7v+j0HsR2kLva61fAl5YVdcPsQ/JvA5BVd1dVbc2ly30iol7qurrw+pDapjZ4dgd+N/A\nHcBm4HDgiKq6c4h9TLT01ogv0iD5K+CFwK1V9ZQ+bf4cOAK4Fziuqi4f9kAlLc28StPDvEqjN8in\niPcA6/ttTHIk8PiqWgu8Cjh1SGOTtHzmVZoe5lUasSUL36q6CFhsCvwlwPuatpcCD0uyz3CGJ2k5\nzKs0PcyrNHrDWDeyH3BD6/ZmYP8h7FfS8JlXaXqYV2nIhrVgev5vP3bq1BjSlDGv0vQwr9IQDeM8\nvjfS+zN8c/Zv7ttOEsMqzVNVoz7Zv3mVdtCk5hXMrLSQhTI7jBnfs+j9WVmSHAJ8vTmtzUID2OHL\nG9/4xp16/DAujsExDLP/MRlJXifhZzQJYxh3/45heGMYk4HzCr7HroYxjLv/1TSGfpac8U3yQXp/\nW3rvJDfQO8fqrk3ITquqc5IcmeRaeidq9sTl0piYV2l6mFdp9JYsfKvq6AHanDic4UjaGeZVmh7m\nVRq9qfnLbTMzM+MegmNwDBPT/zSYhOdo3GMYd/+OYbLGMMkm4flxDOPvvwtjWPIvtw2to6RG1Zc0\nDZJQo/9lmYGYV2l7k5xX2D6zyeiG6euEJlW/zA7jrA6SJGmijKIgndjPAVJfU7PUQZIkSdoZY5nx\n9WsYSZIkjdoYlzr4NYwkSZJGx6UOkiRJ6gQLX0mSJHWCha8kSZI6wcJXkiRJnWDhK0mSpE6w8JUk\nSVInWPhKkiSpEyx8JUmS1AkWvpIkSeoEC19JkiR1goWvJEmSOsHCV5IkSZ1g4StJkqROsPCVJElS\nJ1j4SpIkqRMsfCVJktQJFr6SJEnqBAtfSZIkdYKFryRJkjphycI3yfokm5Jck+SkBbbvneTcJFck\nuTLJcSsyUkkDMbPS9DCv0milqvpvTNYAVwM/CdwIfA44uqo2ttpsAH6gql6XZO+m/T5Vdf+8fdVc\nX0mA/v0OT1js+KRxSkJVZcj7HEpm23mVNNl5bdr5Hiu19MvsUjO+64Brq+q6qtoKnAkcNa/NzcAe\nzfU9gK/ND6SkkTGz0vQwr9KI7bLE9v2AG1q3NwPPmtfm3cD5SW4Cdgf+y/CGJ2mZzKw0PcyrNGJL\nFb6DfIfxeuCKqppJ8jjgY0meVlXfnN9ww4YNrVuzwMyAw5Sm3+zsLLOzsyvdzdAy287rzMwMMzMz\nwxynNNGmLa/ge6y6bdDMLrXG9xBgQ1Wtb26/DthWVW9ttTkH+MOquri5/QngpKq6bN6+XH8ktazQ\nmsGhZNY1vtL2Jjmvzf2+x0otO7rG9zJgbZLHJHkg8DLgrHltNtFbmE+SfYAfBr6880OWtAPMrDQ9\nzKs0Yosudaiq+5OcCJwHrAFOr6qNSU5otp8G/BHwniSfp1dI/3ZV3bHC45a0ADMrTQ/zKo3eoksd\nhtqRX8NI21mJr06HxaUO0vYmOa/ge6w0344udZAkSZJWBQtfSZIkdYKFryRJkjrBwleSJEmdYOEr\nSZKkTrDwlSRJUidY+EqSJKkTLHwlSZLUCRa+kiRJ6gQLX0mSJHWCha8kSZI6wcJXkiRJnWDhK0mS\npE6w8JUkSVInWPhKkiSpEyx8JUmS1AkWvpIkSeoEC19JkiR1goWvJEmSOsHCV5IkSZ1g4StJkqRO\nsPCVJElSJ1j4SpIkqRMsfCVJktQJSxa+SdYn2ZTkmiQn9Wkzk+TyJFcmmR36KCUNzMxqJSUZ2aUL\nzKs0Wqmq/huTNcDVwE8CNwKfA46uqo2tNg8DLgZeUFWbk+xdVbcvsK+a66v3gta/3+EJix2fNE5J\nqKqhvrsPK7PtvEptXX39nuS8Nu18j5Va+mV2qRnfdcC1VXVdVW0FzgSOmtfmGODDVbUZYKFAShoZ\nMytND/MqjdhShe9+wA2t25ub+9rWAg9PckGSy5K8YpgDlLQsZlaaHuZVGrFdltg+yHcYuwI/BhwG\nPAT4dJLPVNU1Ozs4SctmZqXpYV6lEVuq8L0ROKB1+wB6n0jbbgBur6r7gPuSXAg8Dfi+UG7YsKF1\naxaYWd5opSk2OzvL7OzsSncztMy28zozM8PMzMwKDFeaTNOWV/A9Vt02aGaX+uW2XegtvD8MuAn4\nLN+/8P6JwDuBFwA/AFwKvKyqrpq3LxfeSy0r9MsyQ8msv9ymfrr6+j3JeW3a+R4rtfTL7KIzvlV1\nf5ITgfOANcDpVbUxyQnN9tOqalOSc4EvANuAd88PpKTRMLPS9DCv0ugtOuM71I78NCptZyVmkIbF\nGV/109XX70nOK/geK823o6czkyRJklYFC19JkiR1goWvJEmSOsHCV5IkSZ1g4StJkqROsPCVJElS\nJ1j4SpIkqRMsfCVJktQJFr6SJEnqBAtfSZIkdYKFryRJkjrBwleSJEmdYOErSZKkTrDwlSRJUidY\n+EqSJKkTLHwlSZLUCRa+kiRJ6gQLX0mSJHWCha8kSZI6wcJXkiRJnWDhK0mSpE6w8JUkSVInWPhK\nkiSpEyx8JUmS1AkWvpIkSeqEJQvfJOuTbEpyTZKTFmn3n5Lcn+RnhztEScthZqXpYV6l0Vq08E2y\nBngnsB54MnB0kif1afdW4FwgKzBOSQMws9L0MK/S6C0147sOuLaqrquqrcCZwFELtPs14EPAbUMe\nn6TlMbPS9DCv0ogtVfjuB9zQur25ue+7kuxHL6inNnfV0EYnabnMrDQ9zKs0YrsssX2QgL0N+J9V\nVUnCIl/DbNiwoXVrFpgZYPfS6jA7O8vs7OxKdzO0zLbzOjMzw8zMzDDGJ02Facsr+B6rbhs0s6nq\nn7skhwAbqmp9c/t1wLaqemurzZf5XhD3Bu4Ffrmqzpq3r5rrq5fdUXxoDYsdnzROSaiqoa7XG1Zm\n23mV2rr6+j3JeW3a+R4rtfTL7FKF7y7A1cBhwE3AZ4Gjq2pjn/bvAc6uqr9fYJuhlFpW6I10KJm1\n8FU/XX39nuS8Ntt8j5Va+mV20aUOVXV/khOB84A1wOlVtTHJCc3201ZktJJ2iJmVpod5lUZv0Rnf\noXbkp1FpOysxgzQszviqn66+fk9yXsH3WGm+HZrxlaRJ0HsjHx3fzCVpdbLwlTQlRlWMTuykniRp\nJy35J4slSZKk1cDCV5IkSZ1g4StJkqROsPCVJElSJ1j4SpIkqRMsfCVJktQJFr6SJEnqBAtfSZIk\ndYKFryRJkjrBwleSJEmdYOErSZKkTrDwlSRJUidY+EqSJKkTLHwlSZLUCRa+kiRJ6gQLX0mSJHWC\nha8kSZI6wcJXkiRJnWDhK0mSpE6w8JUkSVInWPhKkiSpEyx8JUmS1AkDFb5J1ifZlOSaJCctsP2/\nJvl8ki8kuTjJU4c/VEmDMK/S9DCv0mgtWfgmWQO8E1gPPBk4OsmT5jX7MvATVfVU4GTgXcMeqKSl\nmVdpephXafQGmfFdB1xbVddV1VbgTOCodoOq+nRV3dXcvBTYf7jDlDQg8ypND/Mqjdgghe9+wA2t\n25ub+/o5HjhnZwYlaYeZV2l6mFdpxHYZoE0NurMkzwN+ETh0h0ckaWeYV2l6mFdpxAYpfG8EDmjd\nPoDep9LtNAvu3w2sr6o7F9rRhg0bWrdmgZkBhylNv9nZWWZnZ1e6mxXJ68zMDDMzM8McpzTRpi2v\n4Husum3QzKZq8Q+cSXYBrgYOA24CPgscXVUbW20OBM4Hjq2qz/TZT831lYRlfNDdCWGp45PGJQlV\nlSHvc+h5nQSje80AXzcW19XX70nOa9PO91ippV9ml5zxrar7k5wInAesAU6vqo1JTmi2nwb8HrAX\ncGovcGytqnXDPABJSzOv0vQwr9LoLTnjO7SO/DQqbWclZpCGxRnfyTn2SdPV1+9Jziv4HivN1y+z\n/uU2SZIkdYKFryRJkjrBwleSJEmdYOErSZKkTrDwlSRJUidY+EqSJKkTLHwlSZLUCRa+kiRJ6gQL\nX0mSJHWCha8kSZI6wcJXkiRJnWDhK0mSpE6w8JUkSVInWPhKkiSpEyx8JUmS1AkWvpIkSeoEC19J\nkiR1goWvJEmSOsHCV5IkSZ1g4StJkqROsPCVJElSJ1j4SpIkqRMsfCVJktQJFr6SJEnqhF3GPQCN\nV5KR9VVVI+tLkiRpviVnfJOsT7IpyTVJTurT5s+b7Z9P8ozhDxNgdmV2u5wRzK7WMdQyLxfswGOG\na9w/i3H3v5hJyexkPEez4+19Ap6DSRjDuH8OMCnPw/eblLz6M5qMMYy7/y6MYdHCN8ka4J3AeuDJ\nwNFJnjSvzZHA46tqLfAq4NSVGersyux2OSNY5f8ZBjc77gEM9XlIsuzL8573vB163ErPsE9SZv2/\nOhnPwSSMYdw/B5iU52F7k5RXf0aTMYZx978SYxjle+wglprxXQdcW1XXVdVW4EzgqHltXgK8D6Cq\nLgUelmSf5T0t0rgtdwb7jTvwmJEs9TCz0vQwr6vYjhRub3rTmyZyUmXnjeI9djBLFb77ATe0bm9u\n7luqzf4Dj0DSMJlZaXqY11Vv1UyqrB5V1fcC/Bzw7tbtY4F3zGtzNnBo6/bHgR9bYF878pP04mVV\nXxbL345cGFJmx/28ePEyiZdJzauZ9eJl4ctCuVvqrA43Age0bh9A79PmYm32b+7bTlVN+jy8tBoM\nJbPmVRoJ32OlEVtqqcNlwNokj0nyQOBlwFnz2pwFvBIgySHA16tqy9BHKmkQZlaaHuZVGrFFZ3yr\n6v4kJwLnAWuA06tqY5ITmu2nVdU5SY5Mci1wD/ALKz5qSQsys9L0MK/S6MU/KiBJkqQumIo/WTzI\nCb5XuP+/SrIlyRdH3XdrDAckuSDJl5JcmeTXR9z/g5JcmuSKJFclefMo+583ljVJLk9y9pj6vy7J\nF5oxfHYcY5hk5nX8eW3GMBGZNa+Tr+uZNa/fN5ZVndmJn/FN7wTfVwM/SW9B/+eAo6tq4wjH8Bzg\nbuD9VfWUUfU7bwyPAh5VVVck2Q34F+CnR/w8PKSq7k2yC/Ap4DVV9alR9d8ax28BzwR2r6qXjKH/\nrwDPrKo7Rt33pDOv3x3D2PPajGPsmTWvk83MmtcFxrGqMzsNM76DnOB7RVXVRcCdo+xzgTHcUlVX\nNNfvBjYC+454DPc2Vx9Ibz3ayN9IkuwPHAn8H2Ccv8Xsb1AvzLwyGXlt+h5rZs3rVOh8Zs3r93Qh\ns9NQ+A5ygu9OSfIY4BnApSPu9wFJrgC2ABdU1VWj7L9xCvBaYNsY+p5TwMeTXJbkl8c4jklkXucZ\nV16bvsedWfM6+cxsS8fzCh3I7DQUvpO9FmPEmq9hPgT8RvPJdGSqaltVPZ3eeSR/IsnMKPtP8iLg\n1qq6nPF+Ej20qp4BHAH89+ZrOvWY15Zx5hXGm1nzOjXMbKPLeYXuZHYaCt9BTvDdCUl2BT4MfKCq\n/nFc46iqu4B/Ag4ecdf/GXhJs/7ng8Dzk7x/xGOgqm5u/r0N+Ad6XxWqx7w2JiWvMLbMmtfpYGYx\nr41OZHYaCt9BTvC96iUJcDpwVVW9bQz9753kYc31BwM/BVw+yjFU1eur6oCq+iHg5cD5VfXKUY4h\nyUOS7N5cfyhwODC2swdMIPPK+PPajGGsmTWvU6PzmTWvPV3J7MQXvlV1PzB3gu+rgL8Zw29afhC4\nBHhCkhuSjOME4ofS+zvuz2tO8XF5kvUj7P/RwPnN+qNLgbOr6hMj7H8h4/iKbh/gotbz8JGq+ugY\nxjGRzOt3jTuvMHmZNa8TyMwC5rWfVZnZiT+dmSRJkjQMEz/jK0mSJA2Dha8kSZI6wcJXkiRJnWDh\nK0mSpE6w8JUkSVInWPhKkiSpEyx8p1SS30lyZZLPN+ccXJfkN5oTXy/12P8xSDtJw2NmpelhXlcv\nz+M7hZI8G/gz4LlVtTXJw4EHARcDB1fV15Z4/FcGaSdpOMysND3M6+rmjO90ehRwe1VtBaiqO4Cf\nB/YFLkjyCYAkpyb5XPOpdUNz368v0O7wJJck+Zckf9v8mUCSvCXJl5pPvH8y8qOUVg8zK00P87qK\nOeM7hZrQfAp4CPBxen9i8sLmU+Yzm5CSZK+qujPJmqbdr1XVle12SfYGPgysr6r7kpwEPBD4C+CS\nqnpis689quobIz9YaRUws9L0MK+rmzO+U6iq7gGeCbwKuA34myTHNZvTavqyJP8C/CvwI8CTF9jd\nIc39lyS5HHglcCBwF/DvSU5P8jPAfStxLFIXmFlpepjX1W2XcQ9AO6aqtgGfBD6Z5IvAcXObAJL8\nEPBqeuuM7kryHnprlBbysao6Zv6dSdYBh9H7iufE5rqkHWBmpelhXlcvZ3ynUJInJFnbuusZwHXA\nN4E9mvv2AO4BvpFkH+CIVvt2u0uBQ5M8rtn3Q5Osbb7qeVhV/TPwW8DTVup4pNXOzErTw7yubs74\nTqfdgHckeRhwP3ANva9kjgHOTXJjVR3WfK2yCbiB3nqlOe+a1+444INJfqDZ/jv0gvv/kjyI3lc7\nvzmKA5NWKTMrTQ/zuor5y22SJEnqBJc6SJIkqRMsfCVJktQJFr6SJEnqBAtfSZIkdYKFryRJkjrB\nwleSJEmdYOErSZKkTrDwlSRJUidY+EqSJKkTLHwljUySK5P8xLjHsRxJzknyihXuY0OSM5rrByb5\nZpIMad+nJnlDc30myQ3D2G+zv+ck2TSs/UnSSttl3AOQtHokuRuY+zvoDwX+HfhOc/tVVfWjYxnY\nTqiqI0fRTau/64Hdl3pAkuOA46vqOYvuuOpXd3p03+tzG/D4qvpys++LgCcOa/+StNIsfCUNTVXt\nNnc9yVfoFWbnj3FIO2xuxrWqaqm2kyrJA6pq27B3O+T9SdLIuNRB0sgkuS7J85vrG5L8XZIzknwj\nyReSrE3yuiRbknw1yU+1HrtnktOT3JRkc5KTkyz4GpZkXZLLktyV5JYkf9badkiSS5LcmeSKJM9t\nbZtN8gdJLgbuBh7b3Hd8q80vJrkqyR1Jzk1yYGvbKc3Y72qO50f6jO+HknyyOe6PAnu3tj0myba5\nY0tyXJJ/a9p+OckxSZ4I/CXw7GZZxB1N2/c2SxvOaWbfn9fcd/K8/l+X5LYkX0lyzLzjbx/rcUku\naq5f2Nz9+abPl85fOpHkSc0+7myWtby4te29Sf4iyUeaY/lMkscu97mTpJ1h4StplObPnr4IeD+w\nF3A58LHm/n2Bk4HTWm3fC3wbeBzwDOBw4Jf69PN24JSq2hN4LPC3AEn2Az4C/H5V7QW8Bvhwkke0\nHntss9/dga82Y67m8UcBrwN+hl6xehHwwWbbC4DnAGubfl8KfK3P+P4v8DngEc1x/rcFnhuSPLQ5\nlvVVtQfwbOCKqtoEnAB8uqp2r6qHtx52NHByM/v+qfb4G49q+t236fddSdY22+a3/a6qmlub/dSm\nz7+bN9ZdgbOBc4FHAr8G/HWSJ7SavQzYQO/nfS3wh81jl/PcSdIOs/CVNE4XVtXHquo7wIfoFWRv\naW7/DfCYJHsk2Qc4AvjNqrqvqm4D3ga8vM9+vw2sTbJ3Vd1bVZc29x8LnFNV5wJU1ceBy4AXNtsL\neG9VbayqbVV1/7z9/grw5qq6ullC8Gbg6c2s77fpFctPapYYXF1Vt8wfWNP2YOB3q2prs072bPov\nIdgGPCXJg6tqS1VdNberBdoW8I9V9enm+L7Vp+1c3xcC/0SvIN1ZhwAPraq3VNX9VXUBvQ8ZR7fa\n/H1VXdb8fP8aeHpz/1YGeO4kaWdZ+Eoap1tb1+8Dbm+tqb2v+Xc34CBgV+Dm5mv0O+l91f/IPvs9\nHngCsDHJZ5PMFbYHAS+d20ezn0PpzYLOWeysBwcBb289dm5Wct+m0Hsn8BfAliSnJVnol9T2Be6s\nqvta9311oc6q6h56RemvADc1ywR+eJHxLTV++vT96CUeM4h9F+j7q8390CvKt7S23UfvZ0uzDnyQ\n506SdoqFr6RpcAPwLeARVbVXc9mzqp6yUOOquraqjqmqRwJvBT6U5CHA9cAZrX3s1Xxt/8fthy8y\njuvpnZ2i/fiHVtVnmn7fUVUHA0+mV3i/doF93Azs1YxnzkH9+q2qj1bV4fSK803AuwcY5/ftpnV9\nob5vaq7fQ+9sHHPaHwiWchNwQLLdadgOAm4caICDPXeStFMsfCVNvKq6Gfgo8L+S7J7kAUkelz7n\nBE5ybJK52eC76BV+3wE+ALw4yeFJ1iR5UPMLWvu1H77IUP4SeH2SJzf97Jnkpc31g5M8q1nrei/b\nn8qtfSxfpbe84k1Jdk3y4/TWOi90HD+Y5Khmre9WeoXp3D63APs3/S029ixw/1zfz6G3zGNuve4V\nwM8meXCSx9ObOW/bQm+N9UIupXfcv93se6Y5rjMXGdvccQ703EnSzrLwlTQuC/0i1WK3Xwk8ELgK\nuINesdZvRvIFwJVJvgmcAry8qr5VVZuBo4DX01tmcT3warYvyvrOpFbVP9KbQT4zyV3AF5u+APYA\n3tWM7TrgduBP+uzqGOBZTdvfA97X57gfAPwmvVnTr9H7BbC58/J+AvgScEuSW1uPW+g5bN93M3An\nvRnaM4ATqur/N9tOobdWeQvwHnofFNqP3QC8r1nq8fPtfVfVt4EX01uLfRu9pQuvaO17sZ/3cp47\nSdphmeJTVEqSJEkDc8ZXkiRJnWDhK0mSpE6w8JUkSVInWPhKkiSpEyx8JUmS1AkWvpIkSeoEC19J\nkiR1goWvJEmSOsHCV5IkSZ3wHxnGM7AdTotEAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 16 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot empirical marginal distributions at T=100 with initial states 3, 4, and 5." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Write your own code\n", "inits = [3, 4, 5]\n", "\n", "T = 100\n", "\n", "fig, axes = plt.subplots(1, 3, figsize=(12, 3))\n", "\n", "for init, ax in zip(inits, axes):\n", " draw_histogram(cross_sectional_dist(mc0, init=init, T=T), ax=ax,\n", " title='Initial state = {0}'.format(init), \n", " xlabel='States')\n", "\n", "fig.suptitle('Empirical marginal distribution', y=-0.05, fontsize=12)\n", "plt.show()\n" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAr4AAAD3CAYAAAD7T1f4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH85JREFUeJzt3XuYZHV95/H3xwFUMlwl4U5ICCZqvKDuBMMmtpLVAY2Y\nq8Ggi5LAkyfEZEWX1Zg40UST6K7GkCXETNCgKzGSGDAIRqABAVESMCIDC+osM6BcHBiuRnC++0ed\nlpqmq7t6proufd6v56lnqs751fl9T9V8q77nd351OlWFJEmStNw9YdQBSJIkScNg4StJkqRWsPCV\nJElSK1j4SpIkqRUsfCVJktQKFr6SJElqBQvfIUlyfpLXzLP+9CRv63Nb00lOGFx0krqZr9JkMWfV\nLwvf7ZBkfZIj+2lbVUdX1VnN845Pcvms9b9RVX/YZ9fV3BYT68FJtiTp+z1v9u/Fi+lnEJJckuTO\nJPclWZfk14cdg5Yf83VpJXlhE/M7RxWDlhdzdmk0/T6U5P7mdsGwYxilHUYdwIRbdHKMgSyibS2y\n/aC8Abixqh5Jsgq4LMllVXXTCGLR8mG+LpEkOwJ/BnyeyXuNNb7M2aVRwMur6uIR9D1yjvgOSHOE\n+bkk70myKcnXkqzuWj+d5IQkPwb8JfCC5khrU7P+QzMjJUn2SPKpZtRzU5LzkuzfZxyrklyTZHOS\nbyZ5b7Pqsubfe5t+fyLJIUkuTnJ3kruSfCTJbs12zgIOAs5r2r+pWX54kiuT3JPkuiQvHMTr162q\nvlxVj3QtegC4b9D9qL3M14E7BbgAuIkRFd9a3szZgWttnlr4DtYq4EbgKcCfAmu71hVQVXUjcBJw\nVVXtUlV7dq9v7qd57kHN7WHgtD5j+DPgfVW1G/DDwN83y3+q+Xe3pt+rm8d/BOwLPA04EFhDJ9DX\nALfSOSrcpare23wwfAp4R1XtAbwJOCfJXnMF0nyw3NPjdu58O9E892FgGnh9VX2jz/2X+mW+dtnW\nfE3yg8DrgHfS4i9TDYU522V7vmOBjzaF/4VJntXnvi8LFr6D9f+qam1VFfC3wL5JfmCOdr2+HAJQ\nVZuq6h+r6ttV9QDwLqDfo77vAIcm2auqHupKvsf1WVVfraqLquqRqrobeN8C/RwHnF9VFzTP/yxw\nDXD0XI2r6uVVtUeP2yvm24mqejmwEngt8KEkBy2w39Jima9bb39b8/UDwNuq6kEm89S0Joc5u/X2\ntzVnXw38YHO7BLhwZiS6DSx8B+ubM3eq6qHm7srFbiTJzknOSGcC+mbgUmC3JP2MppwAPBVYl+QL\nSV42Tz97Jzk7ycamn7PoHEn38oPAL3UfVQJHAPv0vXOLUFXfrapPAFcDP7cUfajVzNftlORngZVV\nNTPqFRz11dIxZwegqq6qqv+oqoer6o+Be3lsxHrZ88dto9FrRGRm+Sl0EmtVVd2Z5DnAv9H5Qpl3\nNKWqbqFzNEeSXwA+kWTPHs97F/Bd4Mer6t4krwT+fJ44bwXOqqoT54thRpJPA/+5x+rLqqrnB8Ys\nOwIP9tlWGjTztXe+vhh4fpKZqUi7Ad9N8uNV5cGqRsWcXdx3bKvO0jjiOxp3AAek80voGd0jJSvp\nzDna3CTU2+fYxpxHpkmOS/L9zcPNdP5DbwHuav49pKv5SjoF5X3N3KI3zxFnd/uPAD+b5CVJViR5\nUpKp9PhRQFUd1cxdmus2Z0Im+dEkRyV5cpIdkxwHPB/4zFztpSEwX3t/gf4ecCjwbOA5wLnAX9GZ\n8yuNijnb+zv2wCRHJNmp2f6b6YxCXzFX++XIwndw5prb1uso6iLgK8A3k9w5x/PfDzwZuBu4Evj0\nIrb9UuD6JPfTmU/0K80pjYfoTLK/Ip1fsa4C/gB4Lp3kPQ84Z9Z23w28rTnl8saq2ggcA7wVuJPO\n0ekpDPb/Ueh8CN1B57TWrwEvq6pbB9iHZL4OQFU9UFV3Nrc76BQTD1bVvYPqQ2qYs4OxC/C/gU3A\nRuAlwFFVdc8A+xhr6cwRn6dB8jfAy4A7q+qZPdp8ADgKeAg4vqquHXSgkhZmvkqTw3yVhq+fo4gz\ngdW9ViY5GviRqjoUOBE4fUCxSVo881WaHOarNGQLFr5VdTkw3xD4K4APN22vBnZPsvdgwpO0GOar\nNDnMV2n4BjFvZH9gQ9fjjcABA9iupMEzX6XJYb5KAzaoCdOzf/3YqktjSBPGfJUmh/kqDdAgruN7\nG50/wzfjgGbZVpKYrNIsVTXsi/2br9I2Gtd8BXNWmstcOTuIEd9z6fxZWZIcDtzbXNZmrgB47CoS\ntcjb27fhOY/1OYjb29/+9oFuzxgmN4ZB9D8ii87XSX2PxiGGUfdvDIOLYUT6zlfwO3Y5xDDq/pdT\nDL0sOOKb5GN0/rb0Xkk2NNmxY/Mf/oyqOj/J0UluoXOhZi9cLo2I+SpNDvNVGr4FC9+qOraPNicP\nJhxJ28N8lSaH+SoN3wT95bapUQfA1JQxGMN49D8JxuE1GnUMo+7fGMYrhvE2NeoAxuI9GnUMo+6/\nDTEs+JfbBtZRUjN9JWE4P0zNvPM8pFFKQg3/xzJ96c5XSeOdr+B3rDRbr5ydoBFfSZIkadtZ+EqS\nJKkVLHwlSZLUCha+kiRJagULX0mSJLWCha8kSZJawcJXkiRJrWDhK0mSpFaw8JUkSVIrWPhKkiSp\nFSx8JUmS1AoWvpIkSWoFC19JkiS1goWvJEmSWsHCV5IkSa1g4StJkqRWsPCVJElSK1j4SpIkqRUs\nfCVJktQKFr6SJElqBQtfSZIktYKFryRJklrBwleSJEmtsGDhm2R1khuT3Jzk1DnW75XkgiTXJbk+\nyfFLEqmkvpiz0uQwX6XhSlX1XpmsAG4Cfga4DfgicGxVretqswZ4YlW9JcleTfu9q+rRWduqmb6S\nAL37HZww3/5Jo5SEqsqAtzmQnO3OV0njna9NO79jpS69cnahEd9VwC1Vtb6qHgHOBo6Z1eYbwK7N\n/V2Bb81OSElDY85Kk8N8lYZshwXW7w9s6Hq8EfiJWW0+CFyc5HZgF+CXBxeepEUyZ6XJYb5KQ7ZQ\n4dvPOYy3AtdV1VSSQ4B/SfLsqrp/dsM1a9Z0PZoGpvoMU5p809PTTE9PL3U3A8vZ7nydmppiampq\nkHFKY23S8hX8jlW79ZuzC83xPRxYU1Wrm8dvAbZU1Z90tTkf+KOquqJ5fBFwalVdM2tbzj+SuizR\nnMGB5KxzfKWtjXO+Nsv9jpW6bOsc32uAQ5McnGQn4FXAubPa3EhnYj5J9gZ+FPja9ocsaRuYs9Lk\nMF+lIZt3qkNVPZrkZOBCYAWwtqrWJTmpWX8G8C7gzCRfolNI//eq2rTEcUuagzkrTQ7zVRq+eac6\nDLQjT8NIW1mKU6eD4lQHaWvjnK/gd6w027ZOdZAkSZKWBQtfSZIktYKFryRJklrBwleSJEmtYOEr\nSZKkVrDwlSRJUitY+EqSJKkVLHwlSZLUCha+kiRJagULX0mSJLWCha8kSZJawcJXkiRJrWDhK0mS\npFaw8JUkSVIrWPhKkiSpFSx8JUmS1AoWvpIkSWoFC19JkiS1goWvJEmSWsHCV5IkSa1g4StJkqRW\nsPCVJElSK+ww6gAkSZMjydD6qqqh9SWpHSx8JUmLNIyCdHgFtqT2WHCqQ5LVSW5McnOSU3u0mUpy\nbZLrk0wPPEpJfTNnpclhvkrDlflOJSVZAdwE/AxwG/BF4NiqWtfVZnfgCuClVbUxyV5Vdfcc26qZ\nvjqnyoYzYuCpMo2rJFTVQIe1BpWz3fkqdWvr5/c452vTzu9YqUuvnF1oxHcVcEtVra+qR4CzgWNm\ntXk1cE5VbQSYKyElDY05K00O81UasoUK3/2BDV2PNzbLuh0K7JnkkiTXJHnNIAOUtCjmrDQ5zFdp\nyBb6cVs/5zB2BJ4LHAnsDFyV5PNVdfP2Bidp0cxZaXKYr9KQLVT43gYc2PX4QDpHpN02AHdX1cPA\nw0kuA54NPC4p16xZ0/VoGphaXLTSBJuenmZ6enqpuxlYznbn69TUFFNTU0sQrjSeJi1fwe9YtVu/\nObvQj9t2oDPx/kjgduALPH7i/Y8BpwEvBZ4IXA28qqpumLUtJ95LXZboxzIDyVl/3KZe2vr5Pc75\n2rTzO1bq0itn5x3xrapHk5wMXAisANZW1bokJzXrz6iqG5NcAPw7sAX44OyElDQc5qw0OcxXafjm\nHfEdaEcejUpbWYoRpEFxxFe9tPXze5zzFfyOlWbb1suZSZIkScuCha8kSZJawcJXkiRJrWDhK0mS\npFaw8JUkSVIrWPhKkiSpFSx8JUmS1AoWvpIkSWoFC19JkiS1goWvJEmSWsHCV5IkSa1g4StJkqRW\nsPCVJElSK1j4SpIkqRUsfCVJktQKFr6SJElqBQtfSZIktYKFryRJklrBwleSJEmtYOErSZKkVrDw\nlSRJUitY+EqSJKkVLHwlSZLUCha+kiRJagULX0mSJLXCgoVvktVJbkxyc5JT52n3n5I8muTnBxui\npMUwZ6XJYb5KwzVv4ZtkBXAasBp4OnBskqf1aPcnwAVAliBOSX0wZ6XJYb5Kw7fQiO8q4JaqWl9V\njwBnA8fM0e63gE8Adw04PkmLY85Kk8N8lYZsocJ3f2BD1+ONzbLvSbI/nUQ9vVlUA4tO0mKZs9Lk\nMF+lIdthgfX9JNj7gf9RVZUkzHMaZs2aNV2PpoGpPjYvLQ/T09NMT08vdTcDy9nufJ2ammJqamoQ\n8UkTYdLyFfyOVbv1m7Op6p13SQ4H1lTV6ubxW4AtVfUnXW2+xmOJuBfwEPDrVXXurG3VTF+d3B3G\nQWuYb/+kUUpCVQ10vt6gcrY7X6Vubf38Hud8bdr5HSt16ZWzCxW+OwA3AUcCtwNfAI6tqnU92p8J\nnFdV/zDHOpNS6rJEX6QDyVkLX/XS1s/vcc7XZp3fsVKXXjk771SHqno0ycnAhcAKYG1VrUtyUrP+\njCWJVtI2MWelyWG+SsM374jvQDvyaFTaylKMIA2KI77qpa2f3+Ocr+B3rDRbr5z1L7dJkiSpFSx8\nJUmS1AoLXc5Mkkauc+p2eDx9K0nLk4WvpAkxrGJ0bKdxSpK2k1MdJEmS1AoWvpIkSWoFC19JkiS1\ngoWvJEmSWsHCV5IkSa1g4StJkqRWsPCVJElSK1j4SpIkqRUsfCVJktQKFr6SJElqBQtfSZIktYKF\nryRJklrBwleSJEmtYOErSZKkVrDwlSRJUitY+EqSJKkVLHwlSZLUCha+kiRJagULX0mSJLWCha8k\nSZJaoa/CN8nqJDcmuTnJqXOs/9UkX0ry70muSPKswYcqqR/mqzQ5zFdpuBYsfJOsAE4DVgNPB45N\n8rRZzb4G/HRVPQt4J/BXgw5U0sLMV2lymK/S8PUz4rsKuKWq1lfVI8DZwDHdDarqqqra3Dy8Gjhg\nsGFK6pP5Kk0O81Uasn4K3/2BDV2PNzbLejkBOH97gpK0zcxXaXKYr9KQ7dBHm+p3Y0leBLweOGKb\nI5K0PcxXaXKYr9KQ9VP43gYc2PX4QDpHpVtpJtx/EFhdVffMtaE1a9Z0PZoGpvoMU5p809PTTE9P\nL3U3S5KvU1NTTE1NDTJOaaxNWr6C37Fqt35zNlXzH3Am2QG4CTgSuB34AnBsVa3ranMQcDFwXFV9\nvsd2aqavJCziQHc7hIX2TxqVJFRVBrzNgefrOBjeZwb4uTG/tn5+j3O+Nu38jpW69MrZBUd8q+rR\nJCcDFwIrgLVVtS7JSc36M4DfB/YATu8kHI9U1apB7oCkhZmv0uQwX6XhW3DEd2AdeTQqbWUpRpAG\nxRHf8dn3cdPWz+9xzlfwO1aarVfO+pfbJEmS1AoWvpIkSWoFC19JkiS1goWvJEmSWsHCV5IkSa1g\n4StJkqRWsPCVJElSK1j4SpIkqRUsfCVJktQKFr6SJElqBQtfSZIktYKFryRJklrBwleSJEmtYOEr\nSZKkVrDwlSRJUitY+EqSJKkVLHwlSZLUCha+kiRJagULX0mSJLWCha8kSZJawcJXkiRJrWDhK0mS\npFaw8JUkSVIrWPhKkiSpFSx8JUmS1AoLFr5JVie5McnNSU7t0eYDzfovJTls8GECTC/NZhcTwbQx\nGMN49D+fccnZ8XiNpkfb+xi8BuMQw6jfBxiX1+HxxiVffY/GI4ZR99+GGOYtfJOsAE4DVgNPB45N\n8rRZbY4GfqSqDgVOBE5fmlCnl2azi4lgmf9nMIbJ6b+XccrZ8XiNpkfb+xi8BuMQw6jfBxiX12Fr\n45SvvkfjEcOo+29DDAuN+K4Cbqmq9VX1CHA2cMysNq8APgxQVVcDuyfZe+CRSuqHOStNDvNVGrKF\nCt/9gQ1djzc2yxZqc8D2hyZpG5iz0uQwX6Vhq6qeN+AXgA92PT4O+PNZbc4Djuh6/FnguXNsq7x5\n87b1bb7825YbA8rZUb8u3ryN421c89Wc9eZt7ttcebcD87sNOLDr8YF0jjbna3NAs2wrVZUF+pK0\n/QaSs+arNBR+x0pDttBUh2uAQ5McnGQn4FXAubPanAu8FiDJ4cC9VXXHwCOV1A9zVpoc5qs0ZPOO\n+FbVo0lOBi4EVgBrq2pdkpOa9WdU1flJjk5yC/Ag8Lolj1rSnMxZaXKYr9LwpZkbJEmSJC1rE/GX\n2/q5wPcS9/83Se5I8uVh990Vw4FJLknylSTXJ3nDkPt/UpKrk1yX5IYk7x5m/7NiWZHk2iTnjaj/\n9Un+vYnhC6OIYZyZr6PP1yaGschZ83X8tT1nzdfHxbKsc3bsR3ybC3zfBPwMnQn9XwSOrap1Q4zh\np4AHgL+tqmcOq99ZMewD7FNV1yVZCfwr8Mohvw47V9VDSXYAPge8qao+N6z+u+J4I/A8YJeqesUI\n+v868Lyq2jTsvsed+fq9GEaer00cI89Z83W8mbPm6xxxLOucnYQR334u8L2kqupy4J5h9jlHDN+s\nquua+w8A64D9hhzDQ83dnejMRxv6F0mSA4Cjgb8GRvkrZn9BPTfzlfHI16bvkeas+ToRWp+z5utj\n2pCzk1D49nOB71ZJcjBwGHD1kPt9QpLrgDuAS6rqhmH233gf8GZgywj6nlHAZ5Nck+TXRxjHODJf\nZxlVvjZ9jzpnzdfxZ852aXm+QgtydhIK3/GeizFkzWmYTwC/3RyZDk1Vbamq59C5juRPJ5kaZv9J\nXg7cWVXXMtoj0SOq6jDgKOA3m9N06jBfu4wyX2G0OWu+TgxzttHmfIX25OwkFL79XOC7FZLsCJwD\nfKSqPjmqOKpqM/DPwPOH3PVPAq9o5v98DHhxkr8dcgxU1Teaf+8C/pHOqUJ1mK+NcclXGFnOmq+T\nwZzFfG20ImcnofDt5wLfy16SAGuBG6rq/SPof68kuzf3nwz8F+DaYcZQVW+tqgOr6oeAXwEurqrX\nDjOGJDsn2aW5/33AS4CRXT1gDJmvjD5fmxhGmrPm68Rofc6arx1tydmxL3yr6lFg5gLfNwB/N4Jf\nWn4MuBJ4apINSUZxAfEj6Pwd9xc1l/i4NsnqIfa/L3BxM//oauC8qrpoiP3PZRSn6PYGLu96HT5V\nVZ8ZQRxjyXz9nlHnK4xfzpqvY8icBczXXpZlzo795cwkSZKkQRj7EV9JkiRpECx8JUmS1AoWvpIk\nSWoFC19JkiS1goWvJEmSWsHCV5IkSa1g4TuhkvxukuuTfKm55uCqJL/dXPh6oef+Tj/tJA2OOStN\nDvN1+fI6vhMoyQuA/wm8sKoeSbIn8CTgCuD5VfWtBZ7/9X7aSRoMc1aaHObr8uaI72TaB7i7qh4B\nqKpNwC8C+wGXJLkIIMnpSb7YHLWuaZa9YY52L0lyZZJ/TfLx5s8EkuSPk3ylOeJ9z9D3Ulo+zFlp\ncpivy5gjvhOoSZrPATsDn6XzJyYva44yn9ckKUn2qKp7kqxo2v1WVV3f3S7JXsA5wOqqejjJqcBO\nwF8AV1bVjzXb2rWq7hv6zkrLgDkrTQ7zdXlzxHcCVdWDwPOAE4G7gL9LcnyzOl1NX5XkX4F/A54B\nPH2OzR3eLL8yybXAa4GDgM3At5OsTfJzwMNLsS9SG5iz0uQwX5e3HUYdgLZNVW0BLgUuTfJl4PiZ\nVQBJfgg4hc48o81JzqQzR2ku/1JVr569MMkq4Eg6p3hObu5L2gbmrDQ5zNflyxHfCZTkqUkO7Vp0\nGLAeuB/YtVm2K/AgcF+SvYGjutp3t7saOCLJIc22vy/Joc2pnt2r6tPAG4FnL9X+SMudOStNDvN1\neXPEdzKtBP48ye7Ao8DNdE7JvBq4IMltVXVkc1rlRmADnflKM/5qVrvjgY8leWKz/nfpJO4/JXkS\nnVM7/20YOyYtU+asNDnM12XMH7dJkiSpFZzqIEmSpFaw8JUkSVIrWPhKkiSpFSx8JUmS1AoWvpIk\nSWoFC19JkiS1goWvJEmSWsHCV5IkSa1g4StJkqRWsPCVNFRJ7k9y8DzrT0/ytu3sYyrJhu3ZxlJJ\n8qtJLhzQttYnObLPtscnubzr8bzvwyLjeEuSDzb3D06yJclAvl+SHNTEmkFsT1K77TDqACSNXpL1\nwA8A3+1afGZVvWHQfVXVLgus/41B9zlOquqjwEcHtbnmti1xzPs+QOcAAjirqg5cYFvv3pYYevS5\nHnh9VV3cbPtWYMFYJakfFr6SoFM8vXym2BiVJE+oqi2jjGF7zIxKVtU2FaOTKMmKqvruwi37VoCj\nu5KWhFMdJM2rOUV+RZL/leSeJLck+ckkr0tya5I7kry2q/2Hkvxlks8kuS/JdJKDutZvSfLDXW1P\nT3J+kgeAFzXL3tnV/pgk1yXZ3PT90mb565Lc0PTx1SQnLmKftiT5jSQ3N89/R5JDklyV5N4kZyfZ\nsWm7e5JPJbkzyaYk5yXZv2tb00n+MMkVwIPADyV5SZKbmm39RZJLk5zQ9XpePiuWk5L83+b1Pa1r\n3SFJLk5yd5K7knwkyW597uNTkpzbvG5XA4fM8RrMvA9HJ/lK81psTPLGJDsDnwb2a6Ya3Jdk3yRr\nknwiyVlJNgPHN8vOmhXCCUluS3J7klO6+p39/n5vWkqzjYOA85o+3zR76kSS/Zr9+lbz/v1a17bW\nJPl4kg838V6f5Hn9vF6S2sHCV9KM+UbZVgFfAvYEPgZ8HHgunWLqOOC0plCa8WrgHcBewHXMf2r/\nWOCdVbUS+Bxdp++TrAI+DJxSVbsBPw2sb553B/CyqtoVeB3wviSH9buzwEuAw4DDgVOBDzaxHAQ8\ns7kPnc/Jtc3yg4CHgdNmbes44NeAlcD9wN8329wTuAl4AfNPSXgZ8HzgWcAvzxT3jT8C9gWeBhwI\nrOlz//4CeAjYB3g9ndeoVwxrgROb1/IZwCVV9RCwGri9qnapql2r6htN+1cAf9+8Jx/tsd0p4Efo\nvM6n5rG5yD2nZ1TVa4Bb6Zx92KWq3jtHs7ObNvsCvwi8K8mLutb/LJ3/o7sB5/L490pSi1n4SoJO\n0fvJZsRx5nZC1/qvV9WHm1P4Hwf2A95RVY9U1b8A36FT5Mz4VFV9rqq+A/wu8ILuUdJZPllVVwFU\n1X/MWncCsLaqLmrW315VNzX3z6+qrzf3LwM+A/zUIvb5T6vqgaq6Afgy8OmqWl9V99EZ6Tys2fam\nqvrHqvp2VT0AvAt4Ydd2CvhQVa1rpmkcBVxfVZ+sqi1V9QHgmwvE8sdVdV9VbQAuAZ7T9P3Vqrqo\neZ3vBt43q+85JVkB/Dzw+1X1cFV9hc4BRK+Dm+8Az0iya1VtrqprZzbVo/2VVXVuE+O3e7T7g6bv\n64EzeexAYr7tzivJgcBPAqdW1Xeq6kvAXwOv7Wp2eVVd0Pxf/Qjw7G3pS9LyZOErCTrF2zFVtUfX\nbW3X+ju67j8MUFV3zVq2smtbG7+34aoHgU10iuW5+p3v6gsHAF+da0WSo5J8vjnlfQ9wNPCUebY1\n2+x9mv14ZdPPzknOSOcKCpuBS4Hdkq2uMtC9D/vRtf+N2Y9n6y6MH+rqe+9m2sXGpu+z6G8fv5/O\nbzi647p1nva/QOf1W99M3Th8ge0vtD/M0fdc7/9i7Qdsav5PdW+7+6Cq+318CHhSBnSFCUmTzw8D\nSYMWOqfkOw+SlXRO+d++DdvawNYjyTPbfCJwDvCnwA9U1R7A+SzNj6JOAZ4KrGpO7b+w6ae7r+5T\n97fTKdhnYk334z7NbO9ddK608eNN36+hv8/tu4BH6UzNmHFQj7ZU1TVV9Uo6BfMn6Yzqd8cxO7bZ\ny+dqN7vv25r7DwLd02L26WNbM24H9mz+T3Vvu59CXJIsfCV9zyCLxqOTHJFkJ+CdwFVVddsc7ebq\ns7uoXAu8LsmLkzwhyf5JfhTYqbndDWxJchSduaTbIz3ur6QzArw5yZ7A2xd47j8Dz0znR3k7AL/J\n44u7fuNYSadQvK+ZKvLmfjbQXGXhH4A1SZ6c5OnAf52zs2THdK4tvFvzvPt57LJ2dwBPSbJrj/jm\nW/a2pu9nAMcDf9csv47O/489kuwD/M6s593BrB/ide3XBuBK4N1JnpjkWXTmL39krvaSNJuFr6QZ\nM7+kn7md0yzvd4Sve93/oVMgfovOXNnjejy317YLoKq+SPPDNeBeYBo4qKruB95AZ2RyE535o/+0\nyBjnW9Yd1/uBJ9Mpsq+kM/+35+tRVd8CfonOaPTddH6Udg3wH11tZ/fVK44/oPMjws3AeXRGufu9\nVNrJdArnbwJ/09x69Xsc8PVmOsWJwK82+3IjnR+KfS2dK1rsO0f8vfbpUuAW4LPAe6rqs826s+j8\nUHI9cAGdH6t1P/fddIrme5K8cY5YjwUOpjP6+w905jFf3NVuMf9XJbVMWnS5SUlDkORMYGNV/d6o\nYxkHzfzSDcCrq+rSUccjSW3miK+kQWv9Hx9I5zq+uzdzkd/aLP78KGOSJFn4Shq8bf4zusvIC+ic\n5r+LzjV6XznHpdokSUPmVAdJkiS1giO+kiRJagULX0mSJLWCha8kSZJawcJXkiRJrWDhK0mSpFaw\n8JUkSVIr/H+8YMBeIohbzwAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 17 }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Powers of $P$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The marginal distrubtions at time $T$ are obtained by $P^T$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.set_printoptions(suppress=True) # Suppress printing with floating point notation\n", "\n", "T = 10\n", "print ('P^{0} ='.format(T))\n", "print (np.linalg.matrix_power(mc0.P, T))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "P^10 =\n", "[[ 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0.33398438 0. 0. 0.66601562 0. ]\n", " [ 0.49807228 0.16679179 0.00001694 0.0007677 0.33319974 0.00115155]\n", " [ 0.99902344 0. 0. 0.00097656 0. 0. ]\n", " [ 0. 0.33300781 0. 0. 0.66699219 0. ]\n", " [ 0.99902344 0. 0. 0. 0. 0.00097656]]\n" ] } ], "prompt_number": 14 }, { "cell_type": "code", "collapsed": false, "input": [ "T = 100\n", "print ('P^{0} ='.format(T))\n", "print (np.linalg.matrix_power(mc0.P, T))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "P^100 =\n", "[[ 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0.33333333 0. 0. 0.66666667 0. ]\n", " [ 0.5 0.16666667 0. 0. 0.33333333 0. ]\n", " [ 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0.33333333 0. 0. 0.66666667 0. ]\n", " [ 1. 0. 0. 0. 0. 0. ]]\n" ] } ], "prompt_number": 15 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare the rows with the stationary distributions obtained by `mc0.stationary_distributions`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that `mc0` is aperiodic\n", "(i.e., the least common multiple of the periods of the recurrent class is one),\n", "so that $P^T$ converges as $T \\to \\infty$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "mc0.period" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 16, "text": [ "1" ] } ], "prompt_number": 16 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Exercises 9" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider the Markov chain given by the following stochastic matrix\n", "(where the actual values of non-zero probabilities are not important):" ] }, { "cell_type": "code", "collapsed": false, "input": [ "P = np.zeros((10, 10))\n", "P[0, 3] = 1\n", "P[1, [0, 4]] = 1/2\n", "P[2, 6] = 1\n", "P[3, [1, 2, 7]] = 1/3\n", "P[4, 3] = 1\n", "P[5, 4] = 1\n", "P[6, 3] = 1\n", "P[7, [6, 8]] = 1/2\n", "P[8, 9] = 1\n", "P[9, 5] = 1" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 98 }, { "cell_type": "code", "collapsed": false, "input": [ "np.set_printoptions(precision=3) # Reduce the number of digits printed\n", "print (P)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. ]\n", " [ 0.5 0. 0. 0. 0.5 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. ]\n", " [ 0. 0.333 0.333 0. 0. 0. 0. 0.333 0. 0. ]\n", " [ 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0. 0. 0.5 0. 0.5 0. ]\n", " [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. ]\n", " [ 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. ]]\n" ] } ], "prompt_number": 99 }, { "cell_type": "code", "collapsed": false, "input": [ "mc1 = qe.MarkovChain(P)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 23 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We call the states $0, \\ldots, 9$, respectively, instead of $1, \\ldots, 10$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(a) Check that this Markov chain is irreducible." ] }, { "cell_type": "code", "collapsed": false, "input": [ "mc1.is_irreducible" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 24, "text": [ "True" ] } ], "prompt_number": 24 }, { "cell_type": "markdown", "metadata": {}, "source": [ "(c) Determine the period of this Markov chain." ] }, { "cell_type": "code", "collapsed": false, "input": [ "mc1.period" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 134, "text": [ "3" ] } ], "prompt_number": 134 }, { "cell_type": "markdown", "metadata": {}, "source": [ "(d) Identify the cyclic classes." ] }, { "cell_type": "code", "collapsed": false, "input": [ "mc1.cyclic_classes" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 26, "text": [ "[array([0, 4, 6, 8], dtype=int32),\n", " array([3, 9], dtype=int32),\n", " array([1, 2, 5, 7], dtype=int32)]" ] } ], "prompt_number": 26 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Exercise 11" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us discuss this exercise using the Markov chain from Exercise 9." ] }, { "cell_type": "code", "collapsed": false, "input": [ "print (mc1.P)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. ]\n", " [ 0.5 0. 0. 0. 0.5 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. ]\n", " [ 0. 0.333 0.333 0. 0. 0. 0. 0.333 0. 0. ]\n", " [ 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0. 0. 0.5 0. 0.5 0. ]\n", " [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. ]\n", " [ 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. ]]\n" ] } ], "prompt_number": 96 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Denote the period of $P$ by $d$:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "d = mc1.period" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 135 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Reorder the states so that the transition matrix is in block form:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "permutation = np.concatenate(mc1.cyclic_classes)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 28 }, { "cell_type": "code", "collapsed": false, "input": [ "P = mc1.P[permutation, :][:, permutation]\n", "print (P)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0. 0. 0.333 0.333 0. 0.333]\n", " [ 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. ]\n", " [ 0.5 0.5 0. 0. 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0.5 0.5 0. 0. 0. 0. 0. 0. ]]\n" ] } ], "prompt_number": 136 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us create a new MarkovChain instance:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "mc2 = qe.MarkovChain(P)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 137 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obtain the block components $P_0, \\ldots, P_d$:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "P_blocks = []\n", "\n", "for i in range(d):\n", " P_blocks.append(mc2.P[mc2.cyclic_classes[i%d], :][:, mc2.cyclic_classes[(i+1)%d]])\n", " print ('P_{0} ='.format(i))\n", " print (P_blocks[i])" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "P_0 =\n", "[[ 1. 0.]\n", " [ 1. 0.]\n", " [ 1. 0.]\n", " [ 0. 1.]]\n", "P_1 =\n", "[[ 0.333 0.333 0. 0.333]\n", " [ 0. 0. 1. 0. ]]\n", "P_2 =\n", "[[ 0.5 0.5 0. 0. ]\n", " [ 0. 0. 1. 0. ]\n", " [ 0. 1. 0. 0. ]\n", " [ 0. 0. 0.5 0.5]]\n" ] } ], "prompt_number": 138 }, { "cell_type": "markdown", "metadata": {}, "source": [ "(b) Show that $P^d$ is block diagonal." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hint: You may use [np.linalg.matrix_power](http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.matrix_power.html)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Write your own code\n", "from numpy import linalg as LA\n", "P_power_d=LA.matrix_power(P,d)\n", "print (P_power_d)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 0.167 0.167 0.5 0.167 0. 0. 0. 0. 0. 0. ]\n", " [ 0.167 0.167 0.5 0.167 0. 0. 0. 0. 0. 0. ]\n", " [ 0.167 0.167 0.5 0.167 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.833 0.167 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0. 0. 0.333 0.333 0. 0.333]\n", " [ 0. 0. 0. 0. 0. 0. 0.333 0.333 0. 0.333]\n", " [ 0. 0. 0. 0. 0. 0. 0.333 0.333 0. 0.333]\n", " [ 0. 0. 0. 0. 0. 0. 0.167 0.167 0.5 0.167]]\n" ] } ], "prompt_number": 139 }, { "cell_type": "markdown", "metadata": {}, "source": [ "(c) Verify that the $i$th diagonal block of $P^d$ equals\n", "$P_i P_{i+1} \\ldots P_{d-1} P_0 \\ldots P_{i-1}$.\n", "\n", "Store these diagonal blocks in a list called `P_power_d_blocks`." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Write your own code\n", "P_power_d_0=np.array([P_power_d[0,[0,1,2,3]],P_power_d[1,[0,1,2,3]],P_power_d[2,[0,1,2,3]],P_power_d[3,[0,1,2,3]]])\n", "\n", "P_power_d_1=np.array([P_power_d[4,[4,5]],P_power_d[5,[4,5]]])\n", "\n", "P_power_d_2=np.array([P_power_d[6,[6,7,8,9]],P_power_d[7,[6,7,8,9]],P_power_d[8,[6,7,8,9]],P_power_d[9,[6,7,8,9]]])\n", "\n", "P_power_d_blocks = [P_power_d_0,P_power_d_1,P_power_d_2]\n", "\n", "\n", "for i in range(d):\n", " \n", " print ('P_power_d_{0} ='.format(i))\n", " print (P_power_d_blocks[i])\n" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "P_power_d_0 =\n", "[[ 0.167 0.167 0.5 0.167]\n", " [ 0.167 0.167 0.5 0.167]\n", " [ 0.167 0.167 0.5 0.167]\n", " [ 0. 1. 0. 0. ]]\n", "P_power_d_1 =\n", "[[ 0.833 0.167]\n", " [ 1. 0. ]]\n", "P_power_d_2 =\n", "[[ 0.333 0.333 0. 0.333]\n", " [ 0.333 0.333 0. 0.333]\n", " [ 0.333 0.333 0. 0.333]\n", " [ 0.167 0.167 0.5 0.167]]\n" ] } ], "prompt_number": 140 }, { "cell_type": "markdown", "metadata": {}, "source": [ "(d) Obtain the stationary distributions each associated with the diagonal blocks of $P^d$.\n", "\n", "Store them in a list called `pi_s`." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Write your own code\n", "mc_0 = qe.MarkovChain(P_power_d_0)\n", "mc_1 = qe.MarkovChain(P_power_d_1)\n", "mc_2 = qe.MarkovChain(P_power_d_2)\n", "\n", "pi_s_0 = mc_0.stationary_distributions\n", "pi_s_1 = mc_1.stationary_distributions\n", "pi_s_2 = mc_2.stationary_distributions\n", "\n", "\n", "\n", "pi_s = [pi_s_0,pi_s_1,pi_s_2]\n", "\n", "for i in range(d):\n", " \n", " print ('pi_s_{0} ='.format(i))\n", " print (pi_s[i])\n" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[array([[ 0.143, 0.286, 0.429, 0.143]]), array([[ 0.857, 0.143]]), array([[ 0.286, 0.286, 0.143, 0.286]])]\n", "pi_s_0 =\n", "[[ 0.143 0.286 0.429 0.143]]\n", "pi_s_1 =\n", "[[ 0.857 0.143]]\n", "pi_s_2 =\n", "[[ 0.286 0.286 0.143 0.286]]\n" ] } ], "prompt_number": 144 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Verify that $\\pi^{i+1} = \\pi^i P_i$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Write your own code\n", "np.dot(pi_s_0,P_0)\n", "np.dot(pi_s_1,P_1)\n", "np.dot(pi_s_2,P_2)\n", "\n", "print (np.dot(pi_s_0,P_0))\n", "print (np.dot(pi_s_1,P_1))\n", "print (np.dot(pi_s_2,P_2))\n", "\n" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 0.399 0.131 0.016 0.052]]\n", "[[ 0.536 0. ]]\n", "[[ 0.286 0.286 0.143 0.286]]\n" ] } ], "prompt_number": 143 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obtain the unique stationary distribution $\\pi$ of the original Markov chain." ] }, { "cell_type": "code", "collapsed": false, "input": [ "print (mc2.stationary_distributions[0])" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[ 0.048 0.095 0.143 0.048 0.286 0.048 0.095 0.095 0.048 0.095]\n" ] } ], "prompt_number": 128 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Verify that $\\pi = (\\pi^0, \\ldots, \\pi^d) / d$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Write your own code\n" ], "language": "python", "metadata": {}, "outputs": [ { "ename": "ValueError", "evalue": "could not broadcast input array from shape (4) into shape (1)", "output_type": "pyerr", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mValueError\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[1;31m# Write your own code\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mpi_s\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0md\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[1;31mValueError\u001b[0m: could not broadcast input array from shape (4) into shape (1)" ] } ], "prompt_number": 149 }, { "cell_type": "markdown", "metadata": {}, "source": [ "(e)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Write your answer in a Markdown mode, with providing code if necessary.\n", "\n", "\n" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Powers of $P$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Print ($P^1, P^2, \\ldots, P^d$.)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Write your own code\n", "from numpy import linalg as LA\n", "for i in range(d):\n", " j=i+1\n", " P_power_j=LA.matrix_power(P,j)\n", "\n", " print ('P_power_{0} ='.format(j))\n", " print(P_power_j)\n", "\n", "\n" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "P_power_1 =\n", "[[ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0. 0. 0.333 0.333 0. 0.333]\n", " [ 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. ]\n", " [ 0.5 0.5 0. 0. 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0.5 0.5 0. 0. 0. 0. 0. 0. ]]\n", "P_power_2 =\n", "[[ 0. 0. 0. 0. 0. 0. 0.333 0.333 0. 0.333]\n", " [ 0. 0. 0. 0. 0. 0. 0.333 0.333 0. 0.333]\n", " [ 0. 0. 0. 0. 0. 0. 0.333 0.333 0. 0.333]\n", " [ 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. ]\n", " [ 0.167 0.167 0.5 0.167 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.5 0.5 0. 0. 0. 0. ]]\n", "P_power_3 =\n", "[[ 0.167 0.167 0.5 0.167 0. 0. 0. 0. 0. 0. ]\n", " [ 0.167 0.167 0.5 0.167 0. 0. 0. 0. 0. 0. ]\n", " [ 0.167 0.167 0.5 0.167 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.833 0.167 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0. 0. 0.333 0.333 0. 0.333]\n", " [ 0. 0. 0. 0. 0. 0. 0.333 0.333 0. 0.333]\n", " [ 0. 0. 0. 0. 0. 0. 0.333 0.333 0. 0.333]\n", " [ 0. 0. 0. 0. 0. 0. 0.167 0.167 0.5 0.167]]\n" ] } ], "prompt_number": 170 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Print ($P^{2d}$, $P^{4d}$, and $P^{6d}$.)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Write your own code\n", "from numpy import linalg as LA\n", " \n", "P_power_2d=LA.matrix_power(P,2*d)\n", "P_power_4d=LA.matrix_power(P,4*d)\n", "P_power_6d=LA.matrix_power(P,6*d)\n", " \n", "print ('P_power_2d')\n", "print(P_power_2d)\n", "\n", "print ('P_power_4d')\n", "print(P_power_4d)\n", "\n", "print ('P_power_6d')\n", "print(P_power_6d)\n" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "P_power_2d\n", "[[ 0.139 0.306 0.417 0.139 0. 0. 0. 0. 0. 0. ]\n", " [ 0.139 0.306 0.417 0.139 0. 0. 0. 0. 0. 0. ]\n", " [ 0.139 0.306 0.417 0.139 0. 0. 0. 0. 0. 0. ]\n", " [ 0.167 0.167 0.5 0.167 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.861 0.139 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.833 0.167 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0. 0. 0.278 0.278 0.167 0.278]\n", " [ 0. 0. 0. 0. 0. 0. 0.278 0.278 0.167 0.278]\n", " [ 0. 0. 0. 0. 0. 0. 0.278 0.278 0.167 0.278]\n", " [ 0. 0. 0. 0. 0. 0. 0.306 0.306 0.083 0.306]]\n", "P_power_4d\n", "[[ 0.143 0.286 0.428 0.143 0. 0. 0. 0. 0. 0. ]\n", " [ 0.143 0.286 0.428 0.143 0. 0. 0. 0. 0. 0. ]\n", " [ 0.143 0.286 0.428 0.143 0. 0. 0. 0. 0. 0. ]\n", " [ 0.144 0.282 0.431 0.144 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.857 0.143 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.856 0.144 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0. 0. 0.285 0.285 0.144 0.285]\n", " [ 0. 0. 0. 0. 0. 0. 0.285 0.285 0.144 0.285]\n", " [ 0. 0. 0. 0. 0. 0. 0.285 0.285 0.144 0.285]\n", " [ 0. 0. 0. 0. 0. 0. 0.286 0.286 0.141 0.286]]\n", "P_power_6d\n", "[[ 0.143 0.286 0.429 0.143 0. 0. 0. 0. 0. 0. ]\n", " [ 0.143 0.286 0.429 0.143 0. 0. 0. 0. 0. 0. ]\n", " [ 0.143 0.286 0.429 0.143 0. 0. 0. 0. 0. 0. ]\n", " [ 0.143 0.286 0.429 0.143 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.857 0.143 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.857 0.143 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0. 0. 0.286 0.286 0.143 0.286]\n", " [ 0. 0. 0. 0. 0. 0. 0.286 0.286 0.143 0.286]\n", " [ 0. 0. 0. 0. 0. 0. 0.286 0.286 0.143 0.286]\n", " [ 0. 0. 0. 0. 0. 0. 0.286 0.286 0.143 0.286]]\n" ] } ], "prompt_number": 178 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Print ($P^{10d + 1}, \\ldots, P^{10d + d}$.)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Write your own code\n", "from numpy import linalg as LA\n", "for i in range(d):\n", "\n", " P_power_10di=LA.matrix_power(P,10*d+i)\n", "\n", " print ('P_power_{0} ='.format(10*d+i))\n", " print(P_power_10di)\n" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "P_power_30 =\n", "[[ 0.143 0.286 0.429 0.143 0. 0. 0. 0. 0. 0. ]\n", " [ 0.143 0.286 0.429 0.143 0. 0. 0. 0. 0. 0. ]\n", " [ 0.143 0.286 0.429 0.143 0. 0. 0. 0. 0. 0. ]\n", " [ 0.143 0.286 0.429 0.143 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.857 0.143 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.857 0.143 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0. 0. 0.286 0.286 0.143 0.286]\n", " [ 0. 0. 0. 0. 0. 0. 0.286 0.286 0.143 0.286]\n", " [ 0. 0. 0. 0. 0. 0. 0.286 0.286 0.143 0.286]\n", " [ 0. 0. 0. 0. 0. 0. 0.286 0.286 0.143 0.286]]\n", "P_power_31 =\n", "[[ 0. 0. 0. 0. 0.857 0.143 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.857 0.143 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.857 0.143 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.857 0.143 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0. 0. 0.286 0.286 0.143 0.286]\n", " [ 0. 0. 0. 0. 0. 0. 0.286 0.286 0.143 0.286]\n", " [ 0.143 0.286 0.429 0.143 0. 0. 0. 0. 0. 0. ]\n", " [ 0.143 0.286 0.429 0.143 0. 0. 0. 0. 0. 0. ]\n", " [ 0.143 0.286 0.429 0.143 0. 0. 0. 0. 0. 0. ]\n", " [ 0.143 0.286 0.429 0.143 0. 0. 0. 0. 0. 0. ]]\n", "P_power_32 =\n", "[[ 0. 0. 0. 0. 0. 0. 0.286 0.286 0.143 0.286]\n", " [ 0. 0. 0. 0. 0. 0. 0.286 0.286 0.143 0.286]\n", " [ 0. 0. 0. 0. 0. 0. 0.286 0.286 0.143 0.286]\n", " [ 0. 0. 0. 0. 0. 0. 0.286 0.286 0.143 0.286]\n", " [ 0.143 0.286 0.429 0.143 0. 0. 0. 0. 0. 0. ]\n", " [ 0.143 0.286 0.429 0.143 0. 0. 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.857 0.143 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.857 0.143 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.857 0.143 0. 0. 0. 0. ]\n", " [ 0. 0. 0. 0. 0.857 0.143 0. 0. 0. 0. ]]\n" ] } ], "prompt_number": 183 }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Simulation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us simulate the Markov chain `mc2`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the frequency distribution of visits to the states along a sample path starting at state 0:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "init = 0\n", "draw_histogram(time_series_dist(mc2, init=init),\n", " title='Time series distribution with init={0}'.format(init),\n", " xlabel='States', ylim=(0, 0.35))\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEZCAYAAACD/A7qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG2BJREFUeJzt3Xu4XXV95/H3xwRGuUqRynDRqMR6ecR6eWIcRLfKYOqo\nsdP6IEgplVHaGdSqtQy2lWOdeRQ7jLXVQWwZL4hia5GJHZRL6wFUisQioCQOEaOES0TAIIiamM/8\nsX6HrJzsffY+OZe9z+98Xs+TJ3uv9fut9VtrJ5+99nevtbZsExER9XjEsAcQERGzK8EeEVGZBHtE\nRGUS7BERlUmwR0RUJsEeEVGZBPuIkvQtSS8c9jimQ9Ilkn5njtcxJun88vhxkn4iSbO07HMk/Wl5\n3JF022wstyzvaEnrZ2t5s71+ScskbZc0UCa099Vsto3ZoZzHPhySHgAmdv7ewM+AX5bnb7T9maEM\nbMRJOhM4wvbAbyCSTgZOsX30NPp0gPNtHz7tQTb9t9OM89bd6T/XJG0EXm/7n8vzZcCtwFLb2+dw\nvR1mtl/fCvwxsBfwOeAPbP9i9kZYhxyxD4ntfWzva3tf4PvAKyaeL7RQVzHscczEoEeq013sHCxz\ntpjRHt8uJL0MOB14CfB44InAu4c6qBGVYB9RkjZKekl5PCbp7yWdL+l+STdKWi7pDEmbJX1f0r9v\n9d1f0nmS7pC0SdJ7egWXpBWS1kraIukuSWe35q2U9DVJ90n6pqQXteaNS/pvkr4KPAA8sUw7pdXm\n9ZJulnSvpC9Jelxr3gfK2LeU7Xl6j/E9QdKVZbsvAx7TmrdT+UDSyZK+W9reKukESU8BPgI8v5Rt\n7i1tP15KBJeUT08vLtPeM2n9Z0i6W9L3JJ0wafvb23qypKvL46vK5BvKOl8zubQj6allGfeVstsr\nW/M+LunDkv6xbMu/SHpij/3zCUlvK48PLfvjP5fnT5J0T3n88PrVlLIeB3yhjO+PWos8sfx7ulvS\nO7utszXG97SWvUnS28preoeaT0k7tZW0F/BF4JCy3vslHdxrHV38LvC3ttfZ/jHw58DJU3dZnBLs\no2tyjewVwCeBA4DrgcvL9EOA9wDnttp+HPgF8CTgWcCxwH/qsZ4PAh+wvT/NEdDfQRMSwD8Cf277\nAOCPgH+QdGCr74lluROfOjwxbkmrgTOA36QJ46uBz5R5LwOOBpaX9b4GuKfH+D4NXAccWLbzd7vs\nGyTtXbZlle39gOcD37S9HjgVuKZ8GvqVVrfjgffY3gf4Snv8xcFlvYeU9X5U0vIyb3Lbh9me+G7k\nyLLOv5801j2ALwBfAg4C3gRcIOnJrWbHAWM0r/cG4L933TswDnTK4xfRlFNe2Hp+1eQOpYz1A3Z8\nSvwfrdlHAU8GXgq8q7wxdt1Mdt7+xwL70eyrU4APS9q/3db2T4FVwB1lvfvZvqu8Ad83xZ/DynKe\nBtzQWueNwGMlHdBjjItWgn3huMr25bZ/SVNbPBB4X3n+WWCZpP0kPRb4DeCtth+yfTfwl8Breyz3\nF8BySY+x/VPb15bpJwKX2P4SgO0rgLXAfyjzDXy8HD1tt71t0nJ/H3iv7e+Umu17gV8vR+2/oHkz\neKqkR5Q2d00eWGn7XODPbG+1fTVNIPYqIWwHniHpUbY32755YlFd2hq42PY1Zft+3qPtxLqvAv4v\nTeDO1Epgb9vvs73N9pdp3kSPb7W5yPba8vpeAPx6j2VdBbxAkmjeLN9PE87QBPuV0xzbu23/3PaN\nNCH6zCnatvfVVpqDgF/a/iLNp7hf69J2l9fC9qdtHzDFn02l6T7AllbX+8vf+05j+xaFBPvC8cPW\n44eAH3nHN98Plb/3oak97gHcOXHEQ1OKOKjHck+hOUJbJ+nrkiaC+/HAa9pHTjSB0f7oPNVZI48H\nPtjqO3FEfkgJsg8BHwY2SzpXUrf/nIcA99l+qDXt+91WZvtBmtD9feCOUsb4tW5tBxw/Pdb9b/v0\nGcQhXdb9/TIdmjedza15D9G8truw/V3gQZrgP5rmDeKOcvT/QqYf7O032J/SfLE/iHsmfen6015j\nnoEHaD4VTJj4RPCTWV7Pgpdgr89twM+BA1tHPPvbfka3xrY32D7B9kHAWcDnSi30BzRnL7SPnPa1\n/f529ynG8QOas3va/fe2/S9lvX9t+7k0H6+fDLyjyzLuBA4o45nw+F7rtX2Z7WNp3nzWA38zwDh3\nWUzrcbd131EeP8jOoTedWvEdwOHlKLu97NunsYy2K2nKWXvYvqM8P5mmjPPNHn1m43S43dmv3cpo\nrys1925/7m+VYr7Nzp9cnglstn3f7g2/Xgn2yti+E7gM+J+S9pX0iPIlWtdz4iWdKGniaH4LzX+8\nXwKfAl4p6VhJSyQ9snxJdmi7+xRD+QjwTklPK+vZX9JryuPnSnpeqTX/lJ1P9Wxvy/dpyj/vlrSH\npBfQfNfQbTt+VdLqUmvfShO8E8vcDBxW1jfV2NVl+sS6j6YpQ03Uy78J/EdJj5J0BM0nn7bNNN9x\ndHMtzXb/cVl2p2zXhVOMbSpXAqexo54+Xp5f3fpUN9lU42vrNZZu+2qqZUy03QwcKOnhI2/bF7TO\nCJv8Z79WKeaTwClqvng+APgz4GMDjmFRSbAvDN2+qJvq+UnAnsDNwL00YdTriPJlwLck/QT4APDa\nUmPdBKwG3klTBvoB8HZ2/s/c84jN9sU0nwAulLQFuKmsC5qP0x8tY9sI/Aj4ix6LOgF4Xmn7LuAT\nPbb7EcBbaY5676EpS/xBmfdPNEd7d0n6Yatft33YnnYncB/NEfb5wKm2/1+Z9wGa7wo204TLpyb1\nHQM+UUpRv91edjnv+pU034XcTVOW+p3Wsgd5vduuoil7TAT7V4FHsesXp+1lvBf40zK+t02xjl7r\nnTzGqcbX3vb1NF+i36rmbKmBP+nYvpTmO4Qv0/y7+S5w5qD9F5O+FyhJWkXz5dsSmlONzpo0fzXN\naUfby593tC562EjzBccvga22V8z2BkRExM6mDHZJS4DvAMfQHAldBxxve12rzd7liyskPQP4vO0j\nyvPvAc+xfe/cbUJERLT1K8WsADbY3mh7K00NcHW7wUSoF/vQfKxuW1BXt0VELHT9gv1Qdj4ta1OZ\nthNJr5a0juaqsje3Zhm4Qs2VjW+Y6WAjIqK/pX3mD3Q6U/mi7OJy5sD57Lgw4Sjbd5azLi6XtL5c\nZBIREXOkX7DfDrTvwnY4zVF7V7avlrRU0oG27ymn3mH7bkmfpynt7BTsknJ7yYiI3WC7a6m7Xylm\nLc3l5ssk7UlzZd+adoNyjrTK42eXld0jaa+JqwnLucXH0pzy1m1wM/pz5plnzngZtfzJvsi+yL5Y\nHPtiKlMesdveJuk04FKa0x3Ps71O0qll/rnAbwEnSdpKc8nvxD1JDgYuKpm/FLjA9mVTjiYiImas\nXykGNzf0+eKkaee2Hr+f5qKByf1upfeNiyIiYo5UceVpp9MZ9hBGRvbFDtkXO2Rf7LAY9sXQfxpP\nkoc9hoiIhUYS3s0vTyMiYoFJsEdEVCbBHhFRmQR7RERlEuwREZVJsEdEVCbBHhFRmQR7RERlEuwR\nEZVJsEdEVCbBHhFRmQR7RERlEuwREZVJsEdEVCbBHhFRmQR7RERlEuwREZVJsEdEVCbBHhFRmaXD\nHkBEP1LXn3WcM/kN3ljoEuyxQMxX2M7vm0jEXOhbipG0StJ6SbdIOr3L/NWSbpB0vaRvSHrJoH0j\nImL2aaqPnZKWAN8BjgFuB64Djre9rtVmb9sPlsfPAD5v+4hB+pY+zkffmEpTipm/I/b8e4yFQBK2\nu37E7HfEvgLYYHuj7a3AhcDqdoOJUC/2AX40aN+IiJh9/YL9UOC21vNNZdpOJL1a0jrgi8Cbp9M3\nIiJmV78vTwf6TGr7YuBiSUcD50t6ynQGMTY29vDjTqdDp9OZTveIiOqNj48zPj4+UNt+NfaVwJjt\nVeX5GcB222dN0ee7NGWY5YP0TY09+kmNPWJXM6mxrwWWS1omaU/gOGDNpIU/SeVEY0nPBrB9zyB9\nIyJi9k1ZirG9TdJpwKXAEuA82+sknVrmnwv8FnCSpK3AA8Brp+o7d5sSERHQpxQzLwNIKSb6SCkm\nYlczKcVERMQCk2CPiKhMgj0iojIJ9oiIyiTYIyIqk2CPiKhMgj0iojIJ9oiIyiTYIyIqk2CPiKhM\ngj0iojIJ9oiIyiTYIyIqk2CPiKhMgj0iojIJ9oiIyiTYIyIqk2CPiKhMgj0iojIJ9oiIyiTYIyIq\nk2CPiKhMgj0iojJ9g13SKknrJd0i6fQu818n6QZJN0r6qqQjW/M2lunXS/r6bA8+IiJ2tXSqmZKW\nAB8CjgFuB66TtMb2ulazW4EX2t4iaRXwUWBlmWegY/ve2R96RER00++IfQWwwfZG21uBC4HV7Qa2\nr7G9pTy9Fjhs0jI0KyONiIiB9Av2Q4HbWs83lWm9nAJc0npu4ApJayW9YfeGGBER0zFlKYYmmAci\n6cXA64GjWpOPsn2npIOAyyWtt3315L5jY2MPP+50OnQ6nUFXGxGxKIyPjzM+Pj5QW9m9s1vSSmDM\n9qry/Axgu+2zJrU7ErgIWGV7Q49lnQk8YPvsSdM91RgiJDGNY4yZro38e4yFQBK2u5a6+5Vi1gLL\nJS2TtCdwHLBm0sIfRxPqJ7ZDXdJekvYtj/cGjgVu2v3NiIiIQUxZirG9TdJpwKXAEuA82+sknVrm\nnwu8CzgAOKc5smKr7RXAwcBFZdpS4ALbl83ZlkREBNCnFDMvA0gpJvpIKSZiVzMpxURExAKTYI+I\nqEyCPSKiMgn2iIjKJNgjIiqTYI+IqEyCPSKiMgn2iIjKJNgjIiqTYI+IqEyCPSKiMgn2iIjKJNgj\nIiqTYI+IqEyCPSKiMv1+8zSGqPxIybzIPcgj6pFgH3nzEbjz9wYSEXMvpZiIiMok2CMiKpNgj4io\nTII9IqIyCfaIiMok2CMiKtM32CWtkrRe0i2STu8y/3WSbpB0o6SvSjpy0L4RETH7NNWFKZKWAN8B\njgFuB64Djre9rtXm+cDNtrdIWgWM2V45SN/S37k4prvmAqX5OY99lF+D+dsPMOr7ImKCJGx3vQil\n3xH7CmCD7Y22twIXAqvbDWxfY3tLeXotcNigfSMiYvb1C/ZDgdtazzeVab2cAlyym30jImIW9Lul\nwMCfSSW9GHg9cNR0+46NjT38uNPp0Ol0Bu0aEbEojI+PMz4+PlDbfjX2lTQ181Xl+RnAdttnTWp3\nJHARsMr2hmn2TY29h9TYG6mxR+xqJjX2tcByScsk7QkcB6yZtPDH0YT6iROhPmjfiIiYfVOWYmxv\nk3QacCmwBDjP9jpJp5b55wLvAg4Azim3md1qe0WvvnO4LRERQZ9SzLwMIKWYnlKKaaQUE7GrmZRi\nIiJigUmwR0RUJsEeEVGZBHtERGUS7BERlUmwR0RUJsEeEVGZBHtERGUS7BERlUmwR0RUJsEeEVGZ\nBHtERGUS7BERlUmwR0RUpt9P40VE7KL89sK8ya2UpyfBHhG7af7ukR/Tk1JMRERlEuwREZVJsEdE\nVCbBHhFRmQR7RERlEuwREZVJsEdEVKZvsEtaJWm9pFsknd5l/lMkXSPpZ5LePmneRkk3Srpe0tdn\nc+AREdHdlBcoSVoCfAg4BrgduE7SGtvrWs3uAd4EvLrLIgx0bN87S+ONiIg++h2xrwA22N5oeytw\nIbC63cD23bbXAlt7LCOXjUVEzKN+wX4ocFvr+aYybVAGrpC0VtIbpju4iIiYvn73ipnpzSCOsn2n\npIOAyyWtt3315EZjY2MPP+50OnQ6nRmuNiKiLuPj44yPjw/UVlPdNU3SSmDM9qry/Axgu+2zurQ9\nE3jA9tk9ltV1viTnzm3dNXfQm499o5G+e9787QcY9X0xKvKaDJ8kbHctdfcrxawFlktaJmlP4Dhg\nTa/1TFrpXpL2LY/3Bo4FbprWyCMiYtqmLMXY3ibpNOBSYAlwnu11kk4t88+VdDBwHbAfsF3SW4Cn\nAb8KXFTu27wUuMD2ZXO3KRERAX1KMfMygJRiekopppGP/aMnr8nwzaQUExERC0yCPSKiMgn2iIjK\nJNgjIiqTYI+IqEyCPSKiMgn2iIjKJNgjIiqTYI+IqEyCPSKiMgn2iIjKJNgjIiqTYI+IqEyCPSKi\nMgn2iIjKJNgjIiqTYI+IqEyCPSKiMgn2iIjKJNgjIiqTYI+IqEyCPSKiMgn2iIjK9A12SaskrZd0\ni6TTu8x/iqRrJP1M0tun0zciImafbPeeKS0BvgMcA9wOXAccb3tdq81BwOOBVwP32T570L6lnaca\nw2ImCZiPfSNG+TWYv/0Ao74vRkVek+GThG11m9fviH0FsMH2RttbgQuB1e0Gtu+2vRbYOt2+EREx\n+/oF+6HAba3nm8q0Qcykb0RE7KalfebP5PPPwH3HxsYeftzpdOh0OjNYbUREfcbHxxkfHx+obb8a\n+0pgzPaq8vwMYLvts7q0PRN4oFVjH6hvauy9pcbeSD139OQ1Gb6Z1NjXAsslLZO0J3AcsKbXembQ\nNyIiZsmUpRjb2ySdBlwKLAHOs71O0qll/rmSDqY542U/YLuktwBPs/1At75zuTEREdGnFDMvA0gp\npqeUYhr52D968poM30xKMRERscAk2CMiKpNgj4ioTII9IqIyCfaIiMok2CMiKpNgj4ioTII9IqIy\nCfaIiMok2CMiKpNgj4ioTII9IqIyCfaIiMok2CMiKtPvp/EiYoQ0t8udP7ldbn+j+Jok2CMWnPm7\nD3oMarRek5RiIiIqk2CPiKhMgj0iojIJ9oiIyiTYIyIqk7NiuhjF05ciIgaVYO9ptE5fiogYVN9S\njKRVktZLukXS6T3a/FWZf4OkZ7Wmb5R0o6TrJX19NgceERHdTXnELmkJ8CHgGOB24DpJa2yva7V5\nOXCE7eWSngecA6wssw10bN87J6OPiIhd9DtiXwFssL3R9lbgQmD1pDavAj4BYPta4NGSHtuan1pD\nRMQ86hfshwK3tZ5vKtMGbWPgCklrJb1hJgONiIjB9PvydNBvEHsdlb/A9h2SDgIul7Te9tWTG42N\njT38uNPp0Ol0BlxtRMRiMQ7snJe9aKpT7SStBMZsryrPzwC22z6r1eYjwLjtC8vz9cCLbG+etKwz\ngQdsnz1pukftdL/mdMf5Oyum1/bP3zh6j2EUjMrrMQpGZV+MyjhGwbD2hSRsdz2o7leKWQssl7RM\n0p7AccCaSW3WACeVFa0Efmx7s6S9JO1bpu8NHAvctLubExERg5myFGN7m6TTgEuBJcB5ttdJOrXM\nP9f2JZJeLmkD8CDwe6X7wcBF5WKfpcAFti+bqw2JiIjGlKWYeRlASjEj8HE3H3Vba8u+2LG2Efi3\nOfU4RsFCLMVERMQCM3K3FJjP+7SM8lFARMTuGrlgb8xP+SEiokYpxUREVCbBHhFRmQR7RERlEuwR\nEZVJsEdEVCbBHhFRmQR7RERlEuwREZVJsEdEVCbBHhFRmQR7RERlRvReMTEq5vOmbDDaN2bLvhgt\neT16S7DHAObvXtOjL/titOT16CalmIiIyiTYIyIqk2CPiKhMgj0iojIJ9oiIyiTYIyIq0zfYJa2S\ntF7SLZJO79Hmr8r8GyQ9azp9IyJidk0Z7JKWAB8CVgFPA46X9NRJbV4OHGF7OfBG4JxB+86e8blZ\n7II0PuwBjJDxYQ9ghIwPewAjZHzYA5hz/Y7YVwAbbG+0vRW4EFg9qc2rgE8A2L4WeLSkgwfsO0vG\n52axC9L4sAcwQsaHPYARMj7sAYyQ8WEPYM71C/ZDgdtazzeVaYO0OWSAvhERMcv6Bfug1+surOtt\nIyIq1u9eMbcDh7eeH05z5D1Vm8NKmz0G6At0u5nP7rxPvHvaPaa+idD8vVfN/jhq3Be7O4bsix2y\nL3aocV/s0C/Y1wLLJS0D7gCOA46f1GYNcBpwoaSVwI9tb5Z0zwB9sZ2j/YiIWTRlsNveJuk04FJg\nCXCe7XWSTi3zz7V9iaSXS9oAPAj83lR953JjIiICtJDuMRwREf0t6CtPcwHUDpIOl/RlSd+W9C1J\nbx72mIZJ0hJJ10v6wrDHMkySHi3pc5LWSbq5lEsXJUlnlP8fN0n6tKR/M+wxzZUFG+zzewHUgrAV\neKvtpwMrgf+yyPfHW4Cbmb9fYhhVHwQusf1U4EhgUZZDy3d9bwCebfsZNOXh1w5zTHNpwQY783oB\n1OizfZftb5bHD9D8Bz5kuKMaDkmHAS8H/pZFfCqupP2Bo23/b2i+97K9ZcjDGpb7aQ5+9pK0FNiL\n5oy+Ki3kYB/k4qlFqRydPAu4drgjGZoPAO8Atg97IEP2BOBuSR+T9K+S/kbSXsMe1DDYvhc4G/gB\nzVl6P7Z9xXBHNXcWcrAv9o/YXUnaB/gc8JZy5L6oSHoF8EPb17OIj9aLpcCzgf9l+9k0Z6391+EO\naTgkPQn4Q2AZzSfZfSS9bqiDmkMLOdgHuXhqUZG0B/APwKdsXzzs8QzJvwNeJel7wGeAl0j65JDH\nNCybgE22ryvPP0cT9IvRc4Gv2b7H9jbgIpp/K1VayMH+8MVTkvakuQBqzZDHNDRqLkc7D7jZ9l8O\nezzDYvudtg+3/QSaL8f+2fZJwx7XMNi+C7hN0pPLpGOAbw9xSMO0Hlgp6VHl/8oxNF+uV6nflacj\nKxdA7eIo4ETgRknXl2ln2P7SEMc0ChZ7ye5NwAXl4Oe7lAsIFxvbN5RPbmtpvnv5V+Cjwx3V3MkF\nShERlVnIpZiIiOgiwR4RUZkEe0REZRLsERGVSbBHRFQmwR4RUZkEe1RP0p+UWxnfUG7lu0LSWyQ9\naoC+fzhIu4hRkvPYo2qSnk9z86cX2d4q6VeARwJfBZ5r+54+/b83SLuIUZIj9qjdwcCPyq2dJ+7y\n99s0N4L6sqR/ApB0jqTrypH9WJn25i7tjpX0NUnfkPR3kvYu099XfsThBkl/Me9bGdGSI/aoWgne\nr9Dcf/sK4LO2rypH4s8pQY+kA2zfV37A5QrgTba/1W4n6TE0N1lbZfuh8qtdewIfprnB1FPKsvaz\nff+8b2xEkSP2qJrtB4HnAG8E7gY+K+nkMrt9W9/jJH2D5h4iT6f5Va7JVpbpXyv34zkJeBywBfiZ\npPMk/Sbw0FxsS8SgFuxNwCIGZXs7cCVwpaSbgJMnZgFIegLwdppa+hZJH6Opw3dzue0TJk+UtAJ4\nKU2Z57TyOGIocsQeVZP0ZEnLW5OeBWwEfgLsV6btR/MjFPdLeizwG6327XbXAkeVH21A0t6Slpdy\nz6NtfxF4G/DMudqeiEHkiD1qtw/w15IeDWwDbqEpy5wAfEnS7bZfWkor62l+bvErrf4fndTuZOAz\nrV+4/xOa8P8/kh5JU95563xsWEQv+fI0IqIyKcVERFQmwR4RUZkEe0REZRLsERGVSbBHRFQmwR4R\nUZkEe0REZRLsERGV+f8y0cue8rLncgAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 34 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Observe that the distribution is close to the (unique) stationary distribution." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, plot the simulated marginal distributions\n", "at $T = 10d + 1, \\ldots, 11d, 11d + 1, \\ldots, 12d$ with initial state 0:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "init = 0\n", "T = 10 * d + 1\n", "\n", "fig, axes = plt.subplots(2, d, figsize=(12, 6))\n", "\n", "for t, ax in enumerate(axes.flatten()):\n", " draw_histogram(cross_sectional_dist(mc2, init=init, T=T+t), ax=ax,\n", " title='T = {0}'.format(T+t))\n", "\n", "fig.suptitle('Empirical marginal distributions with init={0}'.format(init), y=0.05, fontsize=12)\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAr4AAAGZCAYAAABvxjJmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X20ZHdd5/v3x26CZCUhwYxRkmAYDA4wREEnRhA4Ckub\n6CWKjtjyMCBK1mgUR8cJQa45PqE4V0UmggFjZBDJOCEXAjeT+ACFgIjkEhKBTiYN5NKdYCAE8gAK\nif29f+zdobpy6lT1OXWqdp39fq21V5+q/Tt7f6tOf8751t6/2pWqQpIkSdruvmrRBUiSJEnzYOMr\nSZKkXrDxlSRJUi/Y+EqSJKkXbHwlSZLUCza+kiRJ6gUbX0mSJPWCje8SSHJ3krva5UCSLw7d3j3D\n/fx2kk8muTPJ/iS/m2Tn0PrXJrk+yb8k+Q+z2q+0nXQhr0kemeStST6d5LNJrkzyyFntW9ouOpLX\nr0ny3iS3JbkjyTVJfmBW+9ahbHyXQFUdVVVHV9XRwP8HfP/B21X1phnu6iLg0VV1DHA68D3ATwyt\n/xDwU8AHAT/5RFpDR/L6YOAtwCOBE4C/B946w31L20JH8no38OPA11bVg4FV4M+THDXD/au1c/IQ\n9UVV3TB0M8AB4FND618NkOSf51yapBHr5bWqPgB84L6VySuBlyU5rqo+N9dCJU3K65eAGwCSfFW7\n7jbgy3Musxc84rsNJXlJks+NWW6f4nvvAvYBb68qjxJJW2hOeX0y8CmbXmlztjKvSa4D/gn4E+AH\nq8rGdwvY+G5DVfVbVXXcmOUhU3zv0cC3As9O8sz5VC3101bnNclJwAXAz2/NI5D6YyvzWlWnAUfT\nTHV4s1MdtoaNr9ZUVdcArwaeu+haJK1vXF6T/CvgL4A/qKr/sYjaJB1qvb+vVfXlqvpvwF3AU+dd\nWx/Y+G5DSV469K7U0eXOw9jUA4AvbFWdkrYur0mOo2l631JVvznruqU+muPf150T1muDbHy3oap6\n+dC7UkeXY9b6njTOTnJs+/XpNFdwuGxozAOSfDXN/5sjknx1ksznUUnb01bkNckxwFXAe6rqpfN7\nNNL2tkV5/fYk35nkiCQPSnIu8NXA383vkfWHja+G/QDwMeAOmkuvvKyqLhta/5fAF4EzgNe2Xz9p\n3kVKAtbP6w8C3wa8YPhoVDvfV9L8rZfXB9LMw78N+CTNm1F3VdXdiyh0u0vV+pdjTfLHwPcBn66q\nx44Z8yrg6TSN0PPb+SuS5sy8SsvDvErzN80R34uBXeNWJjkT+MaqOhV4EfCaGdUm6fCZV2l5mFdp\nziY2vlX1bmC9az8+A3h9O/b9wLFJTphNeZIOh3mVlod5leZvFnN8T6S5GPNB+wHnkUndZF6l5WFe\npRmb1ZvbRt/Zv/7EYUmLZF6l5WFepRnaOYNt3AycPHT7pPa+QyQxrNKIqpr35eDMq7RBXc0rmFlp\nLWtldhZHfC8HngeQ5Azg81V165gCJi7nn3/+VOMWsXS5tq7XZ233XxakN3nten3Wtlz1LcjUeYXJ\nmfXnuj1r63p9XfsbO/GIb5I3AU8Bjk+yDzif5hNHqKoLq+qKJGcm2UvzKSMvOIxQS5oh8yotD/Mq\nzd/Exreqdk8x5pzZlCNpM8yrtDzMqzR/nfvktpWVlUWXMFaXa4Nu12dt21PXn7su12dtG9f1+rqq\n689bl+vrcm3Q7fq6VtvET26b2Y6Smte+pGWQhJr/m2WmYl6lQ3U5r2BmpVHjMtu5I76SJEnSVrDx\nlSRJUi/Y+EqSJKkXbHwlSZLUCza+kiRJ6gUbX0mSJPWCja8kSZJ6wcZXkiRJvWDjK0mSpF6w8ZUk\nSVIv2PhKkiSpF2x8JUmS1As2vpIkSeoFG19JkiT1go2vJEmSesHGV5IkSb1g4ytJkqResPGVJElS\nL9j4SpIkqRd2LroATS/Jpr6/qmZUiSRJ0vKZeMQ3ya4k1ye5Mcm5a6w/PsmVST6U5MNJnr8llapV\nG1zUF2ZWWh7mVZqvrHcUMMkO4AbgacDNwAeA3VW1Z2jMKvDAqjovyfHt+BOq6t6RbZVHHDenOeK7\n0ecwHvHtmCRU1eYO499/mzPJrHmVDtXlvLbjzKw0ZFxmJx3xPR3YW1U3VdU9wCXAWSNjPgUc0359\nDPDZ0UBKmhszKy0P8yrN2aQ5vicC+4Zu7we+fWTM64B3JLkFOBr4kdmVJ+kwmVlpeZhXac4mNb7T\nnDd5KfChqlpJ8gjgL5N8c1XdNTpwdXX1vq9XVlZYWVk5jFKl5TYYDBgMBlu9m5ll1ryqz5Ytr2Bm\n1W/TZnbSHN8zgNWq2tXePg84UFWvGBpzBfAbVfXe9vZfA+dW1dUj23L+0SY5x3d72aI5gzPJrHmV\nDtXlvLb3m1lpyEbn+F4NnJrklCRHAM8CLh8Zcz3NxHySnAB8E/DxzZcsaQPMrLQ8zKs0Z+tOdaiq\ne5OcA1wF7AAuqqo9Sc5u118IvBy4OMm1NI30f6mq27e4bklrMLPS8jCv0vytO9VhpjvyNMymOdVh\ne9mKU6ezYl6lQ3U5r2BmpVEbneogSZIkbQs2vpIkSeoFG19JkiT1go2vJEmSesHGV5IkSb1g4ytJ\nkqResPGVJElSL9j4SpIkqRdsfCVJktQLNr6SJEnqBRtfSZIk9YKNryRJknrBxleSJEm9YOMrSZKk\nXrDxlSRJUi/Y+EqSJKkXbHwlSZLUCza+kiRJ6gUbX0mSJPWCja8kSZJ6wcZXkiRJvTCx8U2yK8n1\nSW5Mcu6YMStJrkny4SSDmVcpaWpmVloe5lWar1TV+JXJDuAG4GnAzcAHgN1VtWdozLHAe4Hvrar9\nSY6vqtvW2Fatty9NlgTY6HMYfP67JQlVlRlvcyaZNa/Sobqc13acmZWGjMvspCO+pwN7q+qmqroH\nuAQ4a2TMjwFvrqr9AGsFUtLcmFlpeZhXac4mNb4nAvuGbu9v7xt2KvCQJO9McnWS586yQEmHxcxK\ny8O8SnO2c8L6ac6bPAB4PPBU4EjgfUn+rqpuHB24urp639crKyusrKxMXai07AaDAYPBYKt3M7PM\nmlf12bLlFcys+m3azE6a43sGsFpVu9rb5wEHquoVQ2POBR5UVavt7T8CrqyqS0e25fyjTXKO7/ay\nRXMGZ5JZ8yodqst5be83s9KQjc7xvRo4NckpSY4AngVcPjLmrcB3JtmR5Ejg24GPzqJoSYfNzErL\nw7xKc7buVIequjfJOcBVwA7goqrak+Tsdv2FVXV9kiuB64ADwOuqylBKC2BmpeVhXqX5W3eqw0x3\n5GmYTXOqw/ayFadOZ8W8Sofqcl7BzEqjNjrVQZIkSdoWbHwlSZLUCza+kiRJ6gUbX0mSJPWCja8k\nSZJ6wcZXkiRJvWDjK0mSpF6w8ZUkSVIv2PhKkiSpF2x8JUmS1As2vpIkSeoFG19JkiT1go2vJEmS\nesHGV5IkSb1g4ytJkqResPGVJElSL9j4SpIkqRdsfCVJktQLNr6SJEnqBRtfSZIk9cLORRcgSZK2\nvySb3kZVzaAS9ZmNryRJmpPNNK6bb5yliVMdkuxKcn2SG5Ocu864f5fk3iTPnG2Jkg6HmZWWh3mV\n5mvdxjfJDuACYBfwaGB3kkeNGfcK4Ep8SSYtjJmVlod5leZv0hHf04G9VXVTVd0DXAKctca4nwEu\nBT4z4/okHR4zKy0P8yrN2aTG90Rg39Dt/e1990lyIk1QX9Pe5cxzaXHMrLQ8zKs0Z5Pe3DZNwF4J\nvKSqKs1bNseehlldXb3v65WVFVZWVqbYvLQ9DAYDBoPBVu9mZpk1r+qzZcsrmFn127SZzXqXBkly\nBrBaVbva2+cBB6rqFUNjPs5Xgng88EXgJ6vq8pFtlZch2Zzmd95Gn8N4GZiOSUJVzXS+3qwya16l\nQ3U5r+24zmd2c3/DwL9jOhzjMjup8d0J3AA8FbgF+Htgd1XtGTP+YuBtVXXZGus6H8qus/HdXrbo\nD+lMMmtepUN1Oa/tus5n1sZX8zQus+tOdaiqe5OcA1wF7AAuqqo9Sc5u11+4JdVK2hAzKy0P8yrN\n37pHfGe6oyV4Ndp1HvHdXrbiCNKsmFfpUF3OKyxHZj3iq3kal9mJH2AhSZIkbQc2vpIkSeoFG19J\nkiT1go2vJEmSesHGV5IkSb1g4ytJkqRemPSRxZIkrau5TNXmeJmq7vHnqu3IxleSNAObuz6rusqf\nq7YXpzpIkiSpF2x8JUmS1As2vpIkSeoFG19JkiT1go2vJEmSesHGV5IkSb1g4ytJkqRe8Dq+ktQz\nfjCBpL6y8ZWkXvKDCST1j1MdJEmS1As2vpIkSeoFG19JkiT1go2vJEmSemGqxjfJriTXJ7kxyblr\nrH92kmuTXJfkvUlOm32pkqZhXqXlYV6l+ZrY+CbZAVwA7AIeDexO8qiRYR8HnlxVpwG/Brx21oVK\nmsy8SsvDvErzN80R39OBvVV1U1XdA1wCnDU8oKreV1V3tDffD5w02zIlTcm8SsvDvEpzNk3jeyKw\nb+j2/va+cV4IXLGZoiRtmHmVhiTZ9LKFzKs0Z9N8gMXUVzlP8l3AjwNPXGv96urqfV+vrKywsrIy\n7aalpTcYDBgMBlu9G/Mq3c9GPqxj0C6/ckgWZmxmeQUzq36b9m9sJn3sZJIzgNWq2tXePg84UFWv\nGBl3GnAZsKuq9q6xnfIjLjenOfKw0ecwfsRoxyShqmZ6OMm8ahqb+10Co79PZr29WZpVbV3Oaztm\n5pnt0/8TbT/jMjvNVIergVOTnJLkCOBZwOUjG38YTSifMy6UkubCvGrpdXx6wiyZV2nOJk51qKp7\nk5wDXAXsAC6qqj1Jzm7XXwj8MnAc8Jr2F849VXX61pUtaS3mVdvH5o4MLgPzKs3fxKkOM9uRp043\nzakO28tWnDqdFfO6vXX9FPYst9flqQ6z5FQH6VDjMjvNm9u0Dc3iVKC/gCRJ0jKx8e217X8qUZIk\n6aCpPrJYkiRJWnY2vpIkSeoFpzpIkqSls9n3qvg+lX6y8ZUkSUtq41c6Uj851UGSJEm9YOMrSZKk\nXrDxlSRJUi/Y+EqSJKkXbHwlSZLUCza+kiRJ6gUvZyZJ0jbhtW2l9dn4SpK0rXhtW2kcpzpIkiSp\nF2x8JUmS1As2vpIkSeoFG19JkiT1go2vJEmSesHGV5IkSb1g4ytJkqRemNj4JtmV5PokNyY5d8yY\nV7Xrr03yuNmXKWlaZlZaHuZVmq91G98kO4ALgF3Ao4HdSR41MuZM4Bur6lTgRcBrNlPQYDDYzLdv\nqS7X1hgsuoCxuvzcdbm2wzXvzHb9uetyfV2urcu/SxqDRRcwE/P/GzvY+LfOxWDRBYzV7bx2u76u\n1TbpiO/pwN6quqmq7gEuAc4aGfMM4PUAVfV+4NgkJ2y0oK49QcO6XFtjsOgCxuryc9fl2jZgrpnt\n+nPX5fq6XFuXf5c0BosuYFbm/Dd2sOFC52Ow6ALG6nZeu11f12qb1PieCOwbur2/vW/SmJM2X5qk\nDTCz0vIwr9KcTWp8p/3A79EP+N7oB4VL2hwzKy0P8yrNW1WNXYAzgCuHbp8HnDsy5g+BHx26fT1w\nwhrbKhcXl0OX9fK3kYUZZXbRz4uLSxeXrubVzLq4rL2slbudrO9q4NQkpwC3AM8Cdo+MuRw4B7gk\nyRnA56vq1tENVdXoK1ZJszeTzJpXaS78GyvN2bqNb1Xdm+Qc4CpgB3BRVe1Jcna7/sKquiLJmUn2\nAl8AXrDlVUtak5mVlod5leYv7SkSSZIkaVvrzCe3TXMR70VJcnKSdyb5SJIPJ/nZRdc0KsmOJNck\neduiaxmW5NgklybZk+Sj7am6zkhyXvtz/Yckf5bkgYuuaVl0NbPmdXO6nFnzunFdzSuY2c3ocl6h\nm5ntROM7zUW8F+we4D9V1WNo3ozw0x2rD+DFwEdpJnR3ye8DV1TVo4DTgD0Lruc+7by6nwQeX1WP\npTnV+KOLrGlZdDyz5nVzOplZ87pxHc8rmNnN6GReobuZ7UTjy3QX8V6YqvrHqvpQ+/XdNP+xHrrY\nqr4iyUnAmcAfcf/L3ixMkgcDT6qqP4ZmPltV3bHgsobdSfML98gkO4EjgZsXW9LS6GxmzevGdTyz\n5nXjOptXMLMb1fG8Qkcz25XGd5qLeHdC+wrmccD7F1vJIX4P+EXgwKILGfFw4DNJLk7ywSSvS3Lk\noos6qKpuB34H+CTNO6o/X1V/tdiqlsZSZNa8HrbOZta8bspS5BXM7GHqbF6hu5ntSuPbtVMHa0py\nFHAp8OL2VenCJfl+4NNVdQ0deiXa2gk8Hnh1VT2e5h3JL1lsSV+R5BHAzwGn0BxdOCrJsxda1PLo\nfGbN64Z0NrPmdVM6n1cwsxvQ2bxCdzPblcb3ZuDkodsn07wi7YwkDwDeDPxpVb1l0fUMeQLwjCSf\nAN4EfHeS/77gmg7aD+yvqg+0ty+lCWlXfBvwt1X12aq6F7iM5vnUZJ3OrHndsC5n1rxuXKfzCmZ2\ng7qcV+hoZrvS+N53Ee8kR9BcxPvyBdd0nyQBLgI+WlWvXHQ9w6rqpVV1clU9nGbS+Duq6nmLrgua\neVvAviSPbO96GvCRBZY06nrgjCQPan/GT6N584Im62xmzevGdTyz5nXjOptXMLMb1fG8QkczO+mT\n2+Zi3EW8F1zWsCcCzwGuS3JNe995VXXlAmsap2untH4GeGP7y/ZjdOji61V1bfvK/WqauVsfBF67\n2KqWQ8cza143p5OZNa8b1/G8gpndjE7mFbqbWT/AQpIkSb3QlakOkiRJ0pay8ZUkSVIv2PhKkiSp\nF2x8JUmS1As2vksgyd1J7mqXA0m+OHR79wz389tJPpnkziT7k/xu+zGDo+Oe19bxwlntW9ouupLX\ndt/DtSz83dRS13QorzuS/HqSm9sxH0zzkcSaMa/qsGTai2i/sKresQXb/ibg5qq6O8lDgb8ALqiq\nPxwacxzwPuDLwCsPfka4pPtbZF6THAAeUVWfmPW+pe1owXn9deAM4AVVtS/Jo4GPVdWXZl1L33Xi\nOr7qhqq6YehmaK6796mRYb8J/D7wI/OqS9L9TZlXz+pJHbBeXtsDSi8GTquqfe34hX/Qw3blL8Vt\nKMlLknxuzHL7FN97F7APeHtVvXVo3ek0H4f4h+O+X9Lh2aq8tv4myaeSvDnJN2zZg5B6Yovy+ljg\nXuDft3m9IclPbfFD6S0b322oqn6rqo4bszxkiu89GvhW4NlJngnN/CPgD4Bzyvkx0sxsRV5bTwa+\nAfg3wC3A29scS9qgLcrrScCDgVOBU4AfBlaTPG0LH0pv2fhqTVV1DfBq4LntXT8FXFdVfz80LHMv\nTNL9rJFXquo9VXVvVd1Bcxr1FJomWNICrZHXf2r//dWq+lJV/QNwCXDmIurb7mx8t6EkLx16V+ro\ncudhbOoBwBfar78b+MH2NMyngCcAv5PkVbOuX+qTLcrr/XYz8q+kDdiivF43ZoxnV7eAje82VFUv\nr6qjxyzHrPU9aZyd5Nj269NpjvJe1g55Ps3Rom8GvgW4GlgFfmnLH5C0jW1FXpM8Osm3tJdIOgr4\nXWA/sGduD0zahrYir1X1MeDdwC8lOSLJo4BnAW+f1+PqExtfDfsB4GPAHcBFwMuq6mAw76iqT7fL\nrTSXM7uzqu5aXLlSr43NK3ACzanSO9oxJwPfX1X/sohCJa2bV4DdNHPyP0vT8L6sqt459yp7YOJ1\nfJP8MfB9wKer6rFjxrwKeDrwReD57fwVSXNmXqXlYV6l+ZvmiO/FwK5xK5OcCXxjVZ0KvAh4zYxq\nk3T4zKu0PMyrNGcTG9+qejfwuXWGPAN4fTv2/cCxSU6YTXmSDod5lZaHeZXmbxZzfE+kuRjzQftp\nrkknqXvMq7Q8zKs0Y7P6yOLRS+Tcb+JwEi/LIY2oqkVcXsq8ShvQ1byCmZXWslZmZ3HE92aadwwf\ndFJ731oFTFzOP//8qcYtYulybV2vz9ruvyxIb/La9fqsbbnqW5Cp8wqTM+vPdXvW1vX6uvY3dhaN\n7+XA8wCSnAF8vprLXUnqHvMqLQ/zKs3YxKkOSd4EPAU4Psk+4HyaTxyhqi6sqiuSnJlkL82nkLxg\nKwuWNJ55lZaHeZXmb2LjW1W7pxhzzmzKgZWVlVltaua6XBt0uz5rmw/zeqgu12dtG9f1+qZlXg/V\n5fq6XBt0u76u1TbxAyxmtqOk5rUvaRkkoRbzZpmJzKt0qC7nFcysNGpcZv3IYkmSJPWCja8kSZJ6\nwcZXkiRJvWDjK0mSpF6w8ZUkSVIv2PhKkiSpF2x8JUmS1As2vpIkSeoFG19JkiT1go2vJEmSesHG\nV5IkSb1g4ytJkqResPGVJElSL9j4SpIkqRdsfCVJktQLNr6SJEnqBRtfSZIk9YKNryRJknrBxleS\nJEm9YOMrSZKkXpjY+CbZleT6JDcmOXeN9ccnuTLJh5J8OMnzt6RSSVMxs9LyMK/SfKWqxq9MdgA3\nAE8DbgY+AOyuqj1DY1aBB1bVeUmOb8efUFX3jmyr1tuX1DdJqKrMeJszyax5lQ7V5by248ysNGRc\nZicd8T0d2FtVN1XVPcAlwFkjYz4FHNN+fQzw2dFAajaSbGpRL5hZaXmYV2nOdk5YfyKwb+j2fuDb\nR8a8DnhHkluAo4EfmV15ur+NvqK38e0JMystD/MqzdmkI77TdFkvBT5UVQ8FvgX4gyRHb7oySRth\nZqXlYV6lOZt0xPdm4OSh2yfTvCId9gTgNwCq6mNJPgF8E3D16MZWV1fv+3plZYWVlZXDLlhaVoPB\ngMFgsNW7mVlmzav6bNnyCmZW/TZtZie9uW0nzUT6pwK3AH/P/Sfe/y5wR1X9SpITgP8XOK2qbh/Z\nlhPvN6mZp7vxqQ4+/92yRW+WmUlmzat0qC7ntR1nZqUh4zK77hHfqro3yTnAVcAO4KKq2pPk7Hb9\nhcDLgYuTXEszdeK/jAZS0nyYWWl5mFdp/tY94jvTHflqdNM84ru9bMURpFkxr9KhupxXMLPSqI1e\nzkySJEnaFmx8JUmS1As2vpIkSeoFG19JkiT1go2vJEmSesHGV5IkSb1g4ytJkqResPGVJElSL9j4\nSpIkqRdsfCVJktQLNr6SJEnqBRtfSZIk9YKNryRJknrBxleSJEm9YOMrSZKkXrDxlSRJUi/Y+EqS\nJKkXbHwlSZLUCza+kiRJ6gUbX0mSJPWCja8kSZJ6wcZXkiRJvTCx8U2yK8n1SW5Mcu6YMStJrkny\n4SSDmVcpaWpmVloe5lWar1TV+JXJDuAG4GnAzcAHgN1VtWdozLHAe4Hvrar9SY6vqtvW2Fatty9N\nlgTY6HMYfP67JQlVlRlvcyaZNa/Sobqc13acmZWGjMvspCO+pwN7q+qmqroHuAQ4a2TMjwFvrqr9\nAGsFUtLcmFlpeZhXac4mNb4nAvuGbu9v7xt2KvCQJO9McnWS586yQEmHxcxKy8O8SnO2c8L6ac6b\nPAB4PPBU4EjgfUn+rqpuHB24urp639crKyusrKxMXai07AaDAYPBYKt3M7PMmlf12bLlFcys+m3a\nzE6a43sGsFpVu9rb5wEHquoVQ2POBR5UVavt7T8CrqyqS0e25fyjTXKO7/ayRXMGZ5JZ8yodqst5\nbe83s9KQjc7xvRo4NckpSY4AngVcPjLmrcB3JtmR5Ejg24GPzqJoSYfNzErLw7xKc7buVIequjfJ\nOcBVwA7goqrak+Tsdv2FVXV9kiuB64ADwOuqylBKC2BmpeVhXqX5W3eqw0x35GmYTXOqw/ayFadO\nZ8W8Sofqcl7BzEqjNjrVQZIkSdoWbHwlSZLUCza+kiRJ6gUbX0mSJPWCja8kSZJ6wcZXkiRJvWDj\nK0mSpF6w8ZUkSVIv2PhKkiSpF2x8JUmS1As2vpIkSeoFG19JkiT1go2vJEmSesHGV5IkSb1g4ytJ\nkqResPGVJElSL9j4SpIkqRdsfCVJktQLNr6SJEnqBRtfSZIk9YKNryRJknphYuObZFeS65PcmOTc\ndcb9uyT3JnnmbEuUdDjMrLQ8+pTXJJtepM1at/FNsgO4ANgFPBrYneRRY8a9ArgS8H+mtCBmVloe\n/cxrbWKRNm/SEd/Tgb1VdVNV3QNcApy1xrifAS4FPjPj+iQdHjMrLQ/zKs3ZpMb3RGDf0O397X33\nSXIiTVBf097lyzJpccystDzMqzRnkxrfaQL2SuAlVVU0p2CW/DSMtNTMrLQ8zKs0ZzsnrL8ZOHno\n9sk0r0iHfStwSTvp/Hjg6UnuqarLRze2urp639crKyusrKwcfsXSkhoMBgwGg63ezcwya17VZ8uW\nVzCz6rdpM5vmReSYlclO4AbgqcAtwN8Du6tqz5jxFwNvq6rL1lhX6+1LkzW/+Db6HAaf/25JQlXN\n9OjNrDJrXqVDdTmv7brOZ3Zzf8PAv2M6HOMyu+4R36q6N8k5wFXADuCiqtqT5Ox2/YVbUq2kDTGz\n0vIwr9L8rXvEd6Y7WoJXo13nEd/tZSuOIM2KeZUO1eW8wnJk1iO+mqdxmfWT2yRJktQLNr6SJEnq\nBRtfSZIk9YKNryRJknrBxleSJEm9MOkDLCRJWlrtBz9silcSkLYPG19J0ja3uUtoSdo+nOogSZKk\nXvCIryRpU5xOIGlZ2PhKkmbA6QSSus+pDpIkSeoFj/hKkqT7cQqLtiMbX0mSNIZTWLS9ONVBkiRJ\nvWDjK0mSpF6w8ZUkSVIv2PhKkiSpF2x8JUmS1As2vpIkSeoFG19JkiT1go2vJEmSesHGV5IkSb0w\nVeObZFeS65PcmOTcNdY/O8m1Sa5L8t4kp82+VEnTMK/S8jCv0nxNbHyT7AAuAHYBjwZ2J3nUyLCP\nA0+uqtOAXwNeO+tCJU1mXqXlYV6l+ZvmiO/pwN6quqmq7gEuAc4aHlBV76uqO9qb7wdOmm2ZkqZk\nXqXlYV6lOZum8T0R2Dd0e3973zgvBK7YTFGSNsy8SsvDvEpztnOKMTXtxpJ8F/DjwBPXWr+6unrf\n1ysrK6ysrEy7aWnpDQYDBoPBVu/GvEozMWiXQ7MwYzPLK5hZ9du0f2NTtX7ukpwBrFbVrvb2ecCB\nqnrFyLgY8e34AAAQx0lEQVTTgMuAXVW1d43t1KR9aX1JOIzfk6Pfjc9/tyShqjLjbZpXTdT8Ltmc\n4f8fm/vdBFv5+2lWtXU5r+2YmWd21j/XLv8/0fYzLrPTTHW4Gjg1ySlJjgCeBVw+svGH0YTyOeNC\nKWkuzKumVJtYNCPmVZqziVMdqureJOcAVwE7gIuqak+Ss9v1FwK/DBwHvKY9knBPVZ2+dWVrs2Z9\nxEfdYF6l5WFeN2ezf8f8G9ZPE6c6zGxHnjrdtFlOdfCU0+JtxanTWTGv21ufTmF3earDLPVxqoPT\n/7SezUx1kCRJkpaeja8kSZJ6wcZXkiRJvTDNdXwlSZob33wraavY+EqSOmhzb4KSpLU41UGSJEm9\nYOMrSZKkXrDxlSRJUi/Y+EqSJKkXbHwlSZLUCza+kiRJ6gUbX0mSJPWC1/GVJGmb2OyHf/jBH9ru\nbHwlSdpWNtq8+sEf2v6c6iBJkqResPGVJElSL9j4SpIkqRdsfCVJktQLNr6SJEnqBRtfSZIk9YKN\nryRJknphYuObZFeS65PcmOTcMWNe1a6/NsnjNlPQYDDYzLdvqS7X1hgsuoCxuvzcdbm2jZhnZrv+\n3HW5vi7X1uXfJY3BoguYmfn+jR1s/FvnYrDoAsbqdl67XV/Xalu38U2yA7gA2AU8Gtid5FEjY84E\nvrGqTgVeBLxmMwV17Qka1uXaGoNFFzBWl5+7Ltd2uOad2a4/d12ur8u1dfl3SWOw6AJmYv5/Ywcb\n/9a5GCy6gLG6nddu19e12iYd8T0d2FtVN1XVPcAlwFkjY54BvB6gqt4PHJvkhJlXKmkaZlZaHuZV\nmrNJje+JwL6h2/vb+yaNOWnzpUnaADMrLQ/zKs1bVY1dgB8CXjd0+znAfxsZ8zbgiUO3/wp4/Brb\nKhcXl0OX9fK3kYUZZXbRz4uLSxeXrubVzLq4rL2slbudrO9m4OSh2yfTvNpcb8xJ7X2HqKpM2Jek\nzZtJZs2rNBf+jZXmbNJUh6uBU5OckuQI4FnA5SNjLgeeB5DkDODzVXXrzCuVNA0zKy0P8yrN2bpH\nfKvq3iTnAFcBO4CLqmpPkrPb9RdW1RVJzkyyF/gC8IItr1rSmsystDzMqzR/aecGSZIkSdtaZz65\nbZqLeC9KkpOTvDPJR5J8OMnPLrqmUUl2JLkmydsWXcuwJMcmuTTJniQfbU/VdUaS89qf6z8k+bMk\nD1x0Tcuiq5k1r5vT5cya143ral7BzG5Gl/MK3cxsJxrfaS7ivWD3AP+pqh4DnAH8dMfqA3gx8FGa\ndzJ2ye8DV1TVo4DTgD0Lruc+SU4BfpLmHdKPpTnV+KOLrGlZdDyz5nVzOplZ87pxHc8rmNnN6GRe\nobuZ7UTjy3QX8V6YqvrHqvpQ+/XdNP+xHrrYqr4iyUnAmcAfAZ15Z2+SBwNPqqo/hmY+W1XdseCy\nht1J8wv3yCQ7gSNZ493SWlNnM2teN67jmTWvG9fZvIKZ3aiO5xU6mtmuNL7TXMS7E9pXMI8D3r/Y\nSg7xe8AvAgcWXciIhwOfSXJxkg8meV2SIxdd1EFVdTvwO8AngVto3i39V4utamksRWbN62HrbGbN\n66YsRV7BzB6mzuYVupvZrjS+XTt1sKYkRwGXAi9uX5UuXJLvBz5dVdfQoVeirZ3A44FXV9Xjad6R\n/JLFlvQVSR4B/BxwCs3RhaOSPHuhRS2PzmfWvG5IZzNrXjel83kFM7sBnc0rdDezXWl8p7mI90Il\neQDwZuBPq+oti65nyBOAZyT5BPAm4LuT/PcF13TQfmB/VX2gvX0pTUi74tuAv62qz1bVvcBlNM+n\nJut0Zs3rhnU5s+Z14zqdVzCzG9TlvEJHM9uVxneai3gvTJIAFwEfrapXLrqeYVX10qo6uaoeTjNp\n/B1V9bxF1wXNvC1gX5JHtnc9DfjIAksadT1wRpIHtT/jp9G8eUGTdTaz5nXjOp5Z87pxnc0rmNmN\n6nheoaOZnfSRxXMx7iLeCy5r2BNpPkP9uiTXtPedV1VXLrCmcbp2SutngDe2v2w/Rocuvl5V17av\n3K+mmbv1QeC1i61qOXQ8s+Z1czqZWfO6cR3PK5jZzehkXqG7mfUDLCRJktQLXZnqIEmSJG0pG19J\nkiT1go2vJEmSesHGV5IkSb1g4ytJkqResPGVJElSL9j4SpIkqRdsfCVJktQLNr6SJEnqBRtfSZIk\n9YKNryRJknrBxleSJEm9YOMrSZKkXrDxlSRJUi/Y+EqSJKkXbHwlSZLUCza+kiRJ6gUbX0mSJPWC\nja8kSZJ6wcZX2qaS3JXklHXWvybJyza5j5Uk+zazja2S5NlJrprRtm5K8tQpxz4/ybuHbq/7czjM\nOs5L8rr261OSHEgyk9/jSR7W1ppZbG+WklyR5LnrrP+TJL825bamfpxdfk4kbYyNrzRHbQP1xfaP\n6cHlVVuxr6o6uqpuWmf9f6yqX9+KfXdBVb2xqr53Vptrl43Use7PAaZ/AVFVv1lVP7mROtbY501J\nvnto259sa93Q49xKVXVmVb0B7v/C4uAQpvz5HM7jHB2bZJDkhYdbf/u9pyR5Z5IvJNkz7QspSbO1\nc9EFSD1TwPdX1TsWWUSSr6qqA4usYTMOHoHrYpO2VZLsqKp/meEmC/BI5uHZzP+3NwHvBXYB3wdc\nmuTUqrptJpVJmopHfKWOaI9kvTfJ7yb5XJK9SZ6Q5AVJPpnk1iTPGxr/J0n+MMlfJLmzPRr1sKH1\nB5L866Gxr2lPGd8NfNfo6eEkZyX5UJI72n1/b3v/C5J8tN3Hx5K86DAe04Ek/zHJje33/2qSRyR5\nX5LPJ7kkyQPasccmeXuSTye5Pcnbkpw4tK1Bkl9P8l7gC8DDk3xPkhvabf1BkncdPCK3xpSDA0nO\nTvK/2+f3gqF1j0jyjiS3JflMkj9N8uApH+PXJLm8fd7eDzxijefg4M/hzCQfaZ+L/Ul+PsmRwP8C\nHtqeAbgzydcnWU1yaZI3JLkDeH573xtGSnhhkpuT3JLkF4b2O/rzve+ocruNhwFva/f5nzMydSLJ\nQ9vH9dn25/cTQ9taTfLnSV7f1vvhJN86tP7c9vHdmeT6DB1ZHhrz8CSfG7r9uiS3Dt1+Q5IXt18P\nkrwwyb8B/hD4jrbu24c2+ZD2/8+dSf7u4HO+xn5HH+eg/X/5nvZ7r0ryNSNjdyT5DeBJwAU5zDM1\nSR4JPA44v6q+VFWXAdcBPzTtNiTNho2vNH/rHWU7HbgWeAjNEaI/Bx5P00w9h+aP7pFD438M+FXg\neOBDwBvX2fZu4Neq6ijgPQydHk5yOvB64Beq6sHAk4Gb2u+7Ffi+qjoGeAHwe0keN+2DBb6H5o/+\nGcC5wOvaWh4GPLb9GprfRxe19z8M+CfggpFtPQf4CeAo4C7gf7bbfAhwA/AdrH9U7vuAbwNOA34k\nbXPf+g3g64FHAScDq1M+vj8Avgh8HfDjNM/RuBouAl7UPpePAd5ZVV+kOQp4S3ta/Ziq+lQ7/hnA\n/2x/Jm8cs90V4Btpnudz85VT6GNP/1fVc4FP0px9OLqq/q81hl3Sjvl64IeBlyf5rqH1/wfN/9EH\nA5fT/qySfBPw08C3tY/ze/jK/6XhGj4B3Dn0f+nJwF1tc3vw9mD4sVTV9cDZwPvauh/Srg/wozQ/\ns+OAvTQ/z2ntBp4PfC1wBPCf719u/RLwbuCn233/bPt4r2tfSK21HPz/+xjg41X1haFtXtveL2mO\nbHyl+QrwlpE/jsNzBj9RVa9vT+H/OfBQ4Fer6p6q+kvgyzRNzkFvr6r3VNWXgV+iORJ2Imt7S1W9\nD6CqvjSy7oXARVX11+36W6rqhvbrK9omhar6G+AvaI58Teu3q+ruqvoo8A/A/6qqm6rqTpojnY9r\nt317Vf3fVfXPVXU38HLgKUPbKeBPqmpPO03j6cCHq+otVXWgql4F/OOEWn6rqu6sqn3AO4Fvaff9\nsar66/Z5vg34vZF9rynJDuCZwC9X1T9V1UdoXkCMe3HzZeAxSY6pqjuq6pqDmxoz/m+r6vK2xn8e\nM+5X2n1/GLiYr7yQWG+760pyMvAE4Nyq+nJVXQv8EfC8oWHvrqor2/+rfwp8c3v/vwAPpHmcD2jn\nyX58zK7eBawk+Tqan++lwFOSPBw4pt3v/cpb474CLquqq9vpIG+k/dlOoYCLq2pv+xz/+YTvPWT/\nVXVaVR03ZjmnHXYUcMfIdu4Ejp6yRkkzYuMrzVcBZ438cbxoaP2tQ1//E0BVfWbkvqOGtrX/vg03\nR5Nup2mW19rvem+eOgn42Forkjy9PXX82fbU9JnA16yzrVGjj2n09lHtfo5McmGaN13dQdMUPTg5\n5B31w4/hoQw9/tbo7VHDjfEXh/Z9QpppF/vbfb+B6R7jv6J5r8RwXZ9cZ/wP0Tx/N7Wn2M+YsP1J\nj4c19r3Wz/9wPRS4feQI5SeB4RdVwz/HLwJfnWbu+F7g52iOvt6a5E1Jvn7Mft5Fc8T6ScDftLef\nQnO0d/QNbJOs+f9qSsP/LyZ970bm+d4NHDNy37E0za+kObLxlZZXaE7JNzeSo2hO+d+ygW3t49Aj\nyQe3+UDgzcBvA19bVccBV7A1b4r6BeCRwOntqf2ntPsZ3tdw03ELTcN+sNYM357Swe29nOZI5b9t\n9/1cpvv9+BngXpqpGQc9bMxY2iOSP0DTML+F5ujicB2jtY3ev9a40X3f3H79BWB4WszXTbGtg26h\nmTM73AA+jOkacarqTVX1JOAb2v28YszQd9E0vSs00xreAzyR5mc/GLf5aWrYIvfbdztn+64xy6vb\nYR8B/vXI8/nN7f2S5sjGV5q/WTaNZyZ5YpIjgF+jmft48xrj1trncFN5EfCCJN+d5KuSnNjO1Tyi\nXW4DDiR5Os2czc3ImK+PojnadkeShwDnT/je/wd4bJo35e2kmVc62txNW8dRNI3ine1UkV+cZgPt\nafXLgNUkD0ryaOA/rLmz5AFpri384Pb77qJptqE5Wvk1SYaPCo77mY16Wbvvx9DMU/0f7f0fovn/\ncVw7leDnRr7vVkbeiDf0uPYBfwv8ZpIHJjmNZv7yn641fuRxPrL9f/RA4EvAPw89ztH97G3XPwd4\nV1XdBXya5sj4u8bs4lbgpLRvijy420l1TSp7ynH3e86q6jHtnN+1lp9qx/xvmp/H+Um+OskzgX9L\n86JS0hzZ+Erzd/Cd9AeXg3/8pj3CN7zuz2gaxM/SzJV9zpjvHbftAqiqD9C+cQ34PM3Rtoe1jcjP\n0hyZvJ1m/uhbD7PG9e4bruuVwINomuy/pZn/O/b5qKrPAv+e5mj0bTRvSruaptka3fZatQzf/hWa\nNxHeAbyNpiGZ9sjiOTSN8z8Cf9wu4/b7HOAT7XSKFwHPbh/L9TRvFPt4mitafP0a9Y97TO+ieTPX\nXwH/tar+ql33Bpo3UN0EXEnzZrXh7/1Nmqb5c0l+fo1adwOn0Bz9vYxmHvM7hsaNez4f2G77M8Cn\naN54eR7jDYDbhl6wDdp/Pzhm/F/THCn9xySfnqKetaw3dr3/N78P/HD7M3rlOttfy4/SvLHydpo3\n3v1Q+39Y0hyl+nMZTGlbSXIxsL+q/s9F19IF7eWp9gE/VlXjjhZKknrMI77S8ur9hw+kuY7vse1p\n9Ze2d//dImuSJHWXja+0vDb8MbrbyHfQnOb/DM01en9gjUu1SZIEONVBkiRJPeERX0mSJPWCja8k\nSZJ6wcZXkiRJvWDjK0mSpF6w8ZUkSVIv2PhKkiSpF/5/1ZoJF0jysUoAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 35 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare the rows of $P^{10d + 1}, \\ldots, P^{10d + d}$." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }