week12/Lab12.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "# Read files, Linear algebra, and  Markov Chain simulation\n",
    "##### Keywords: Linear algebra, Markov Chain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "## store current working directory\n",
    "cwd      = os.getcwd() \n",
    "\n",
    "## where we store data files\n",
    "datapath = cwd + u'/data/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now to access data stored in the subfolder data\n",
    "## You cannot access houston_weather.txt simply by :\n",
    "## np.loadtxt('weather_comma.txt', delimiter=',')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## You cannot access houston_weather.txt simply by \n",
    "## np.loadtxt('houston_weather.txt', delimiter=',')\n",
    "weather = np.loadtxt(datapath + 'weather_space.txt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear algebra and Markov Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.5   0.25  0.25]\n",
      " [ 0.5   0.    0.5 ]\n",
      " [ 0.25  0.25  0.5 ]]\n"
     ]
    }
   ],
   "source": [
    "weather = np.loadtxt(datapath + 'weather_space.txt')\n",
    "print(weather)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The following are NOT matrix square or matrix product\n",
    "## They are just entrywise operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This is NOT a matrix product, but entrywise squared!!!\n",
      " [[ 0.25    0.0625  0.0625]\n",
      " [ 0.25    0.      0.25  ]\n",
      " [ 0.0625  0.0625  0.25  ]]\n",
      "\n",
      "This is NOT a matrix product, but entrywise product!!!\n",
      " [[ 0.25    0.0625  0.0625]\n",
      " [ 0.25    0.      0.25  ]\n",
      " [ 0.0625  0.0625  0.25  ]]\n",
      "\n",
      "This is NOT a matrix inverse, but entrywise inverse!!!\n",
      " [[  2.   4.   4.]\n",
      " [  2.  inf   2.]\n",
      " [  4.   4.   2.]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/ipykernel_launcher.py:3: RuntimeWarning: divide by zero encountered in reciprocal\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n"
     ]
    }
   ],
   "source": [
    "print(\"This is NOT a matrix product, but entrywise squared!!!\\n\", weather**2)\n",
    "print(\"\\nThis is NOT a matrix product, but entrywise product!!!\\n\", weather*weather)\n",
    "print(\"\\nThis is NOT a matrix inverse, but entrywise inverse!!!\\n\", weather**(-1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# A Matrix is a 2D array that stores real or complex numbers. \n",
    "# You can use numpy to create matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A1 = np.array(    \n",
    "    [[0, 4],\n",
    "     [2, 0]]\n",
    ")\n",
    "\n",
    "A2 = np.array(\n",
    "    [[-1, 2,2],\n",
    "     [1, -2,2]]\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 2) (2, 3)\n"
     ]
    }
   ],
   "source": [
    "print(A1.shape, A2.shape)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## As matrices, A1*A2 is defined, and the product should be a 2 by 3 matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 4, -8,  8],\n",
       "        [-2,  4,  4]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## A1*A2 this will not Run!!\n",
    "\n",
    "matrix1 = np.matrix(\n",
    "    [[0, 4],\n",
    "     [2, 0]]\n",
    ")\n",
    "matrix2 = np.matrix(\n",
    "    [[-1, 2,2],\n",
    "     [1, -2,2]]\n",
    ")\n",
    "\n",
    "matrix1*matrix2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## To take product between matrices A and B\n",
    "1. where A is of shape (m by n) and B is of shape (k by h)\n",
    "1. A*B is defined if and only if n=k\n",
    "l. matrix2 $*$ matrix1 does not run is because you cannot take product between them\n",
    "1. as matrix2 is of shape (2 by 3) and matrix1 is of shape (2 by 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "shapes (2,3) and (2,2) not aligned: 3 (dim 1) != 2 (dim 0)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-13-1049e2f03172>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m## This doesn't run because the shapes of matrices do not match\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mmatrix2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mmatrix1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/matrixlib/defmatrix.py\u001b[0m in \u001b[0;36m__mul__\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m    307\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndarray\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtuple\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    308\u001b[0m             \u001b[0;31m# This promotes 1-D vectors to row vectors\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 309\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mN\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0masmatrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    310\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0misscalar\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mhasattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'__rmul__'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    311\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mN\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: shapes (2,3) and (2,2) not aligned: 3 (dim 1) != 2 (dim 0)"
     ]
    }
   ],
   "source": [
    "## This doesn't run because the shapes of matrices do not match\n",
    "matrix2*matrix1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# if one of the \"arrays\" is of matrix type, you are doing matrix operation!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 4 -8  8]\n",
      " [-2  4  4]]\n",
      "[[ 4 -8  8]\n",
      " [-2  4  4]]\n"
     ]
    }
   ],
   "source": [
    "print(matrix1 * matrix2)\n",
    "print(matrix1 * A2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# np.linalg -- linalg is the linear algebra modulo in python\n",
    "# linalg: https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.linalg.html\n",
    "\n",
    "1. from np import linalg as lin\n",
    "1. lin.det -- determinant\n",
    "1. eigenvalues\n",
    "1. etc.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 4]\n",
      " [2 0]]\n",
      "determinant(matrix1) -8.0\n"
     ]
    }
   ],
   "source": [
    "from numpy import linalg as lin\n",
    "print(matrix1)\n",
    "det = lin.det(matrix1)\n",
    "print('determinant(matrix1)', det)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Eigenvalues & Eigenvectors\n",
    "Let A be an n x n matrix. The number λ is an eigenvalue of A if there exists a non-zero vector C such that\n",
    "\n",
    "Av = λv\n",
    "\n",
    "In this case, vector v is called an eigenvector of A corresponding to λ. You can use numpy to calculate the eigenvalues and eigenvectors of a matrix:\n",
    "\n",
    "The input matrix can either be matrix type or 2d array type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 4]\n",
      " [2 0]]\n",
      "[ 2.82842712 -2.82842712]\n",
      "[ 2.82842712 -2.82842712]\n"
     ]
    }
   ],
   "source": [
    "print(A1)\n",
    "eigvals_a = lin.eigvals(A1)\n",
    "print( eigvals_a)\n",
    "\n",
    "eigvals_m = lin.eigvals(matrix1)\n",
    "print(eigvals_m)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 8.831+0.j    -0.000+0.j     0.584+0.581j  0.584-0.581j]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([\n",
    " [1,2,3,4],   \n",
    " [0,2,3,4],  \n",
    " [0,0,3,4], \n",
    " [1,2,3,4],\n",
    "])\n",
    "\n",
    "eigvals_a = np.round(lin.eigvals(A),3)\n",
    "print(eigvals_a)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve linear equations Ax = b\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.5 ,  0.25])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = matrix1\n",
    "b = np.array([1,1])\n",
    "solu = lin.solve(A,b)\n",
    "solu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## If A is invertable, then x is $A^{-1}b$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 0.5 ],\n",
       "        [ 0.25]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "invA = A**(-1)\n",
    "\n",
    "## But you cannot take invA*b\n",
    "## you can see that the solu above and invA*b here have different shapes of out put\n",
    "## even if the values agree.\n",
    "\n",
    "b = b.reshape((2,1))\n",
    "invA*b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now we want to test what we have seen about markov chain step matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "weather = np.loadtxt(datapath + 'weather_comma.txt',delimiter=',')\n",
    "P = np.matrix(weather)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "P to the power  2\n",
      "[[ 0.4375  0.1875  0.375 ]\n",
      " [ 0.375   0.25    0.375 ]\n",
      " [ 0.375   0.1875  0.4375]]\n",
      "if equals previous power, False\n",
      "\n",
      "P to the power  3\n",
      "[[ 0.40625   0.203125  0.390625]\n",
      " [ 0.40625   0.1875    0.40625 ]\n",
      " [ 0.390625  0.203125  0.40625 ]]\n",
      "if equals previous power, False\n",
      "\n",
      "P to the power  4\n",
      "[[ 0.40234375  0.19921875  0.3984375 ]\n",
      " [ 0.3984375   0.203125    0.3984375 ]\n",
      " [ 0.3984375   0.19921875  0.40234375]]\n",
      "if equals previous power, False\n",
      "\n",
      "P to the power  5\n",
      "[[ 0.40039062  0.20019531  0.39941406]\n",
      " [ 0.40039062  0.19921875  0.40039062]\n",
      " [ 0.39941406  0.20019531  0.40039062]]\n",
      "if equals previous power, False\n",
      "\n",
      "P to the power  6\n",
      "[[ 0.40014648  0.19995117  0.39990234]\n",
      " [ 0.39990234  0.20019531  0.39990234]\n",
      " [ 0.39990234  0.19995117  0.40014648]]\n",
      "if equals previous power, False\n",
      "\n",
      "P to the power  7\n",
      "[[ 0.40002441  0.20001221  0.39996338]\n",
      " [ 0.40002441  0.19995117  0.40002441]\n",
      " [ 0.39996338  0.20001221  0.40002441]]\n",
      "if equals previous power, False\n",
      "\n",
      "P to the power  8\n",
      "[[ 0.40000916  0.19999695  0.3999939 ]\n",
      " [ 0.3999939   0.20001221  0.3999939 ]\n",
      " [ 0.3999939   0.19999695  0.40000916]]\n",
      "if equals previous power, False\n",
      "\n",
      "P to the power  9\n",
      "[[ 0.40000153  0.20000076  0.39999771]\n",
      " [ 0.40000153  0.19999695  0.40000153]\n",
      " [ 0.39999771  0.20000076  0.40000153]]\n",
      "if equals previous power, False\n"
     ]
    }
   ],
   "source": [
    "power_list = list()\n",
    "Pk = np.eye(3)\n",
    "for power in np.arange(0,11):\n",
    "    Pk = P*Pk## If you don't do rounding, you will not see powers agree\n",
    "    power_list.append(Pk)\n",
    "\n",
    "for power in np.arange(2,10):\n",
    "    print(\"\\nP to the power \", power)\n",
    "    print(power_list[power-1])\n",
    "    print('if equals previous power,',\n",
    "          np.array_equal(power_list[power-1] , power_list[power-2]  )\n",
    "         )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "P to the power  2\n",
      "[[ 0.44  0.19  0.38]\n",
      " [ 0.38  0.25  0.38]\n",
      " [ 0.38  0.19  0.44]]\n",
      "if equals previous power, False\n",
      "\n",
      "P to the power  3\n",
      "[[ 0.41  0.2   0.4 ]\n",
      " [ 0.41  0.19  0.41]\n",
      " [ 0.4   0.2   0.41]]\n",
      "if equals previous power, False\n",
      "\n",
      "P to the power  4\n",
      "[[ 0.41  0.2   0.4 ]\n",
      " [ 0.4   0.2   0.4 ]\n",
      " [ 0.4   0.2   0.41]]\n",
      "if equals previous power, False\n",
      "\n",
      "P to the power  5\n",
      "[[ 0.4  0.2  0.4]\n",
      " [ 0.4  0.2  0.4]\n",
      " [ 0.4  0.2  0.4]]\n",
      "if equals previous power, False\n",
      "\n",
      "P to the power  6\n",
      "[[ 0.4  0.2  0.4]\n",
      " [ 0.4  0.2  0.4]\n",
      " [ 0.4  0.2  0.4]]\n",
      "if equals previous power, True\n",
      "\n",
      "P to the power  7\n",
      "[[ 0.4  0.2  0.4]\n",
      " [ 0.4  0.2  0.4]\n",
      " [ 0.4  0.2  0.4]]\n",
      "if equals previous power, True\n",
      "\n",
      "P to the power  8\n",
      "[[ 0.4  0.2  0.4]\n",
      " [ 0.4  0.2  0.4]\n",
      " [ 0.4  0.2  0.4]]\n",
      "if equals previous power, True\n",
      "\n",
      "P to the power  9\n",
      "[[ 0.4  0.2  0.4]\n",
      " [ 0.4  0.2  0.4]\n",
      " [ 0.4  0.2  0.4]]\n",
      "if equals previous power, True\n"
     ]
    }
   ],
   "source": [
    "power_list = list()\n",
    "Pk = np.eye(3)\n",
    "for power in np.arange(0,11):\n",
    "    Pk = np.round(P*Pk,2) ## added rounding and then you can see powers agree\n",
    "    power_list.append(Pk)\n",
    "\n",
    "for power in np.arange(2,10):\n",
    "    print(\"\\nP to the power \", power)\n",
    "    print(power_list[power-1])\n",
    "    print('if equals previous power,',\n",
    "          np.array_equal(power_list[power-1] , power_list[power-2]  )\n",
    "         )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Now let's work with the rat's moving "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "r = np.matrix(pd.read_table(datapath+'rat.txt',delimiter=',',header = None).values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Consider a rat in a maze with 3 cells, indexed as 1,2,3 and the outside (freedom), indexed by 4 (that can only be reached via cell 3). \n",
    "### The rat starts initially in a given cell and then takes a move to another cell, continuing to do so until finally reaching freedom. \n",
    "\n",
    "## Goal: simulate rat moving around and count steps it takes to get to freedom\n",
    "<img src=\"4rat.png\" alt=\"Drawing\" style=\"width: 300px;\"/>\n",
    "\n",
    "$$\n",
    "P = \\begin{pmatrix}\n",
    "0 & 1 & 0 & 0\\\\\n",
    "1/3 & 0 & 1/3 & 1/3\\\\\n",
    "0 & 1 & 0 & 0\\\\\n",
    "0 & 0 & 0 & 1\\\\\n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "S = [1,2,3,4]\n",
    "import random\n",
    "import numpy as np\n",
    "## start moving from location i, \n",
    "## number of steps it takes to get to freedom\n",
    "## recursion\n",
    "def steps(i):\n",
    "    #print(i,'->')\n",
    "    if(i==1):\n",
    "        next_move = 2\n",
    "        return 1+steps(next_move)\n",
    "    \n",
    "    if(i == 2):\n",
    "        next_move = random.choice([1,3,4])\n",
    "        return 1+steps(next_move)\n",
    "        \n",
    "    if(i == 3):\n",
    "        next_move = 2\n",
    "        return 1+steps(next_move)\n",
    "    \n",
    "    if(i == 4):\n",
    "        return 1\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "start at cell 2\n",
      "2 ->\n",
      "3 ->\n",
      "2 ->\n",
      "4 ->\n",
      "4\n",
      "start at cell 2\n",
      "2 ->\n",
      "4 ->\n",
      "2\n",
      "start at cell 2\n",
      "2 ->\n",
      "1 ->\n",
      "2 ->\n",
      "4 ->\n",
      "4\n",
      "start at cell 2\n",
      "2 ->\n",
      "3 ->\n",
      "2 ->\n",
      "3 ->\n",
      "2 ->\n",
      "4 ->\n",
      "6\n",
      "start at cell 2\n",
      "2 ->\n",
      "1 ->\n",
      "2 ->\n",
      "4 ->\n",
      "4\n",
      "start at cell 2\n",
      "2 ->\n",
      "1 ->\n",
      "2 ->\n",
      "4 ->\n",
      "4\n",
      "start at cell 2\n",
      "2 ->\n",
      "4 ->\n",
      "2\n",
      "start at cell 2\n",
      "2 ->\n",
      "4 ->\n",
      "2\n",
      "start at cell 2\n",
      "2 ->\n",
      "1 ->\n",
      "2 ->\n",
      "4 ->\n",
      "4\n",
      "start at cell 2\n",
      "2 ->\n",
      "1 ->\n",
      "2 ->\n",
      "1 ->\n",
      "2 ->\n",
      "1 ->\n",
      "2 ->\n",
      "1 ->\n",
      "2 ->\n",
      "1 ->\n",
      "2 ->\n",
      "4 ->\n",
      "12\n"
     ]
    }
   ],
   "source": [
    "for run in range(10):\n",
    "    print('start at cell',2)\n",
    "    print(steps(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## I want to count in general, \n",
    "## on average how many steps it takes to start from cell 2 and ends with freedom\n",
    "<img src=\"4rat.png\" alt=\"Drawing\" style=\"width: 300px;\"/>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Run N = 100 times, each time compute average of 50 moves\n",
    "\n",
    "def steps(i): \n",
    "    #print(i,'->') turn off print to run simulations\n",
    "    if(i==1):\n",
    "        next_move = 2\n",
    "        return 1+steps(next_move)\n",
    "    \n",
    "    if(i == 2):\n",
    "        next_move = random.choice([1,3,4])\n",
    "        return 1+steps(next_move)\n",
    "        \n",
    "    if(i == 3):\n",
    "        next_move = 2\n",
    "        return 1+steps(next_move)\n",
    "    \n",
    "    if(i == 4):\n",
    "        return 0\n",
    "    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "start move from cell 1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADahJREFUeJzt3X2MZfVdx/H3R5ZqS0mh7mRdge3wByHBxAJOCJWGoNs2\nPBhoY2MgsUIfsqigRZuYtX9oY/9Zkz4YH0LdFiwqYCsFuxZaIdikaaIbF7rCwtqAdGnBhV1KBKrG\nuvD1jznUYZzZe2bu08yP9yuZzL3nnNnz3V+G99w5e+8lVYUkaf37oWkPIEkaDYMuSY0w6JLUCIMu\nSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUiA2TPNnGjRtrdnZ2kqeUpHXvvvvue6aqZgYdN9Ggz87O\nsmfPnkmeUpLWvSSP9znOSy6S1AiDLkmNMOiS1AiDLkmNMOiS1IiBQU9ySpKvJnk4yUNJPtht/0iS\nJ5Ps7T4uHv+4kqTl9Hna4hHgQ1V1f5LjgfuS3NPt+2RVfWx840mS+hoY9Ko6CBzsbr+QZD9w0rgH\nkyStzIquoSeZBc4Cdnebrk3yQJIbk5w44tkkSSvQ+5WiSV4PfAG4rqqeT3I98FGgus8fB963xNdt\nA7YBbNmyZRQza8xmt985tXMf2HHJ1M4trXe9HqEnOZb5mN9cVbcDVNXTVfViVb0EfBo4Z6mvraqd\nVTVXVXMzMwPfikCStEp9nuUS4AZgf1V9YsH2zQsOexewb/TjSZL66nPJ5TzgPcCDSfZ22z4MXJHk\nTOYvuRwArh7LhJKkXvo8y+XrQJbYddfox5EkrZavFJWkRhh0SWqEQZekRhh0SWqEQZekRhh0SWqE\nQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZek\nRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0\nSWrEwKAnOSXJV5M8nOShJB/str8xyT1JHuk+nzj+cSVJy+nzCP0I8KGqOgM4F7gmyRnAduDeqjoN\nuLe7L0makoFBr6qDVXV/d/sFYD9wEnAZcFN32E3AO8c1pCRpsBVdQ08yC5wF7AY2VdXBbtdTwKaR\nTiZJWpHeQU/yeuALwHVV9fzCfVVVQC3zdduS7Emy5/Dhw0MNK0laXq+gJzmW+ZjfXFW3d5ufTrK5\n278ZOLTU11bVzqqaq6q5mZmZUcwsSVpCn2e5BLgB2F9Vn1iwaxdwZXf7SuCLox9PktTXhh7HnAe8\nB3gwyd5u24eBHcDnk7wfeBz4hfGMKEnqY2DQq+rrQJbZvXW040iSVstXikpSIwy6JDXCoEtSIwy6\nJDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDWiz7stSs2b3X7n1M59YMclUzu3\n2uIjdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYY\ndElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYMDHqSG5McSrJvwbaPJHkyyd7u4+Lx\njilJGqTPI/TPAhcusf2TVXVm93HXaMeSJK3UwKBX1deAZycwiyRpCMNcQ782yQPdJZkTRzaRJGlV\nNqzy664HPgpU9/njwPuWOjDJNmAbwJYtW1Z5ulen2e13TnuEiXs1/p2lUVnVI/SqerqqXqyql4BP\nA+cc5didVTVXVXMzMzOrnVOSNMCqgp5k84K77wL2LXesJGkyBl5ySXIrcAGwMckTwO8CFyQ5k/lL\nLgeAq8c4oySph4FBr6orlth8wxhmkSQNwVeKSlIjDLokNcKgS1IjVvs89FcVnxstaT3wEbokNcKg\nS1IjDLokNcKgS1IjDLokNcKgS1IjDLokNcKgS1IjDLokNcKgS1IjDLokNcKgS1IjDLokNcKgS1Ij\nDLokNcKgS1IjDLokNcKgS1IjDLokNcKgS1IjDLokNcKgS1IjDLokNcKgS1IjDLokNcKgS1IjDLok\nNcKgS1IjBgY9yY1JDiXZt2DbG5Pck+SR7vOJ4x1TkjRIn0fonwUuXLRtO3BvVZ0G3NvdlyRN0cCg\nV9XXgGcXbb4MuKm7fRPwzhHPJUlaodVeQ99UVQe7208Bm0Y0jyRplYb+R9GqKqCW259kW5I9SfYc\nPnx42NNJkpax2qA/nWQzQPf50HIHVtXOqpqrqrmZmZlVnk6SNMhqg74LuLK7fSXwxdGMI0larT5P\nW7wV+Afg9CRPJHk/sAN4e5JHgLd19yVJU7Rh0AFVdcUyu7aOeBZJ0hB8pagkNcKgS1IjDLokNcKg\nS1IjDLokNcKgS1IjDLokNcKgS1IjDLokNcKgS1IjDLokNcKgS1IjDLokNcKgS1IjDLokNcKgS1Ij\nDLokNcKgS1IjDLokNcKgS1IjDLokNcKgS1IjDLokNcKgS1IjDLokNcKgS1IjDLokNcKgS1IjDLok\nNcKgS1IjDLokNcKgS1IjDLokNWLDMF+c5ADwAvAicKSq5kYxlCRp5YYKeudnquqZEfw5kqQheMlF\nkhox7CP0Au5OUsCfVtXOxQck2QZsA9iyZcuQp5PaM7v9zmmPMHEHdlwy7RGaNOwj9LdW1dnARcA1\nSc5ffEBV7ayquaqam5mZGfJ0kqTlDBX0qnqy+3wIuAM4ZxRDSZJWbtVBT3JckuNfvg28A9g3qsEk\nSSszzDX0TcAdSV7+c26pqq+MZCpJ0oqtOuhV9Rjw5hHOIkkagk9blKRGGHRJaoRBl6RGGHRJaoRB\nl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RGGHRJaoRBl6RG\nGHRJaoRBl6RGDPM/iZ6o2e13TnsESVrTfIQuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w\n6JLUiHXzwiJJGtY0X6B4YMclYz+Hj9AlqREGXZIaYdAlqRFDBT3JhUm+meTRJNtHNZQkaeVWHfQk\nxwB/AlwEnAFckeSMUQ0mSVqZYR6hnwM8WlWPVdX3gb8CLhvNWJKklRom6CcB31lw/4lumyRpCsb+\nPPQk24Bt3d3vJfnmmE61EXhmTH/2ajlTf2txLmfqZ8Uz5ffHNMkrram16v7Oq53pTX0OGiboTwKn\nLLh/crftFapqJ7BziPP0kmRPVc2N+zwr4Uz9rcW5nKmftTgTrM25xj3TMJdc/gk4LcmpSV4DXA7s\nGs1YkqSVWvUj9Ko6kuRa4O+AY4Abq+qhkU0mSVqRoa6hV9VdwF0jmmVYY7+sswrO1N9anMuZ+lmL\nM8HanGusM6WqxvnnS5ImxJf+S1Ij1l3QkxyT5BtJvrTEvquSHE6yt/v4wIRmOpDkwe6ce5bYnyR/\n2L1FwgNJzl4DM12Q5LkFa/U7E5jphCS3JfmXJPuTvGXR/omvU8+5JrpWSU5fcK69SZ5Pct2iYya6\nVj1nmsb31G8keSjJviS3JvmRRft/OMnnunXanWR23DP1nGs8raqqdfUB/CZwC/ClJfZdBfzxFGY6\nAGw8yv6LgS8DAc4Fdq+BmS5Yag3HPNNNwAe6268BTpj2OvWca+JrteDcxwBPAW9aC2s1YKaJrhPz\nL2T8FvDa7v7ngasWHfOrwKe625cDn1sjc42lVevqEXqSk4FLgM9Me5YVugz485r3j8AJSTZPe6hJ\nSvIG4HzgBoCq+n5V/fuiwya+Tj3nmqatwL9W1eOLtk/ze2q5maZhA/DaJBuA1wH/tmj/Zcz/wAa4\nDdiaJGtgrrFYV0EH/gD4LeCloxzz892voLclOeUox41SAXcnua97Zexi03ibhEEzAbwlyT8n+XKS\nnxjzPKcCh4E/6y6ZfSbJcYuOmcY69ZkLJrtWC10O3LrE9mm+9cZyM8EE16mqngQ+BnwbOAg8V1V3\nLzrsB+tUVUeA54AfXQNzwRhatW6CnuTngENVdd9RDvtbYLaqfhK4h//7yTxub62qs5l/58lrkpw/\nofMezaCZ7mf+V+Y3A38E/M2Y59kAnA1cX1VnAf8BrIW3XO4z16TXCoDMv2DvUuCvJ3G+PgbMNNF1\nSnIi84/ATwV+HDguyS+O85x99JxrLK1aN0EHzgMuTXKA+Xd2/Nkkf7nwgKr6blX9d3f3M8BPTWKw\n7icyVXUIuIP5d6JcqNfbJExypqp6vqq+192+Czg2ycYxjvQE8ERV7e7u38Z8SBea+Dr1mWsKa/Wy\ni4D7q+rpJfZNY62OOtMU1ultwLeq6nBV/Q9wO/DTi475wTp1lz/eAHx3jDP1mmtcrVo3Qa+q366q\nk6tqlvlf+f6+ql7xU2/RNcRLgf3jnivJcUmOf/k28A5g36LDdgG/1D0z4VzmfwU7OM2ZkvzYy9cS\nk5zD/PfC2L7Rq+op4DtJTu82bQUeXnTYRNep71yTXqsFrmD5SxsTX6tBM01hnb4NnJvkdd15t/L/\n/5vfBVzZ3X43890Y94tvBs41rlaN/d0Wxy3J7wF7qmoX8OtJLgWOAM8y/y/J47YJuKP7Pt4A3FJV\nX0nyywBV9SnmX017MfAo8J/Ae9fATO8GfiXJEeC/gMsn8I3+a8DN3a/tjwHvnfI69Z1r4mvV/SB+\nO3D1gm1TXaseM010napqd5LbmL/UcwT4BrBzURNuAP4iyaPMN+Hycc2zwrnG0ipfKSpJjVg3l1wk\nSUdn0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEf8LVEioE2g9pjUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10652cf60>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.0048\n",
      "start move from cell 2\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADWBJREFUeJzt3X2MZfVdx/H3pyw+USzUnawbHjpNJU3Q2AUnaxsaQott\noDTQRqKQWGmD2WpAadLErP1Dq3/RxLY+xGC2gF2V1ja02LVQLEGSSmLRWVx52jYQsk0hCzsU5UFN\nm4Wvf8xZHMeZvffOvXMffrxfyc2995zfnfPJL7ufOXPmnDOpKiRJs+81kw4gSRoNC12SGmGhS1Ij\nLHRJaoSFLkmNsNAlqREWuiQ1wkKXpEZY6JLUiC3j3NjWrVtrfn5+nJuUpJm3f//+Z6pqrte4sRb6\n/Pw8i4uL49ykJM28JN/pZ5yHXCSpERa6JDXCQpekRljoktQIC12SGmGhS1IjLHRJaoSFLkmNsNAl\nqRFjvVJUmlbzu2+f2LYPXX/JxLattriHLkmNsNAlqREWuiQ1wkKXpEZY6JLUCAtdkhphoUtSIyx0\nSWqEhS5JjbDQJakRPQs9yRlJ7knySJKHk1zXLf94kieTHOge79n8uJKk9fRzL5ejwEer6v4kJwP7\nk9zVrft0Vf3h5sWTJPWrZ6FX1WHgcPf6hSQHgdM2O5gkaTADHUNPMg+cA9zXLbo2yQNJbk5y6oiz\nSZIG0HehJ3kt8CXgI1X1PHAD8CZgB8t78J9c53O7kiwmWVxaWhpBZEnSWvoq9CQnslzmt1TVlwGq\n6umqeqmqXgY+A+xc67NVtaeqFqpqYW5ublS5JUmr9HOWS4CbgINV9akVy7evGPZ+4KHRx5Mk9auf\ns1zOAz4APJjkQLfsY8CVSXYABRwCPrwpCSVJfennLJd7gayx6o7Rx5EkbZRXikpSIyx0SWqEhS5J\njbDQJakRFrokNcJCl6RGWOiS1AgLXZIaYaFLUiMsdElqhIUuSY2w0CWpERa6JDXCQpekRljoktQI\nC12SGmGhS1IjLHRJakQ/f1NUGpv53bdPOoI0s9xDl6RGWOiS1AgLXZIaYaFLUiMsdElqhIUuSY2w\n0CWpERa6JDXCQpekRljoktQIC12SGtGz0JOckeSeJI8keTjJdd3y1ye5K8mj3fOpmx9XkrSefvbQ\njwIfraqzgbcC1yQ5G9gN3F1VZwF3d+8lSRPSs9Cr6nBV3d+9fgE4CJwGXAbs7YbtBd63WSElSb0N\ndAw9yTxwDnAfsK2qDnerngK2rfOZXUkWkywuLS0NEVWSdDx9F3qS1wJfAj5SVc+vXFdVBdRan6uq\nPVW1UFULc3NzQ4WVJK2vr0JPciLLZX5LVX25W/x0ku3d+u3Akc2JKEnqRz9nuQS4CThYVZ9asWof\ncFX3+irgK6OPJ0nqVz9/gu484APAg0kOdMs+BlwPfDHJ1cB3gF/anIiSpH70LPSquhfIOqsvHG0c\nSdJGeaWoJDXCQpekRljoktQIC12SGmGhS1IjLHRJaoSFLkmNsNAlqREWuiQ1wkKXpEZY6JLUCAtd\nkhphoUtSIyx0SWqEhS5JjbDQJakRFrokNcJCl6RGWOiS1AgLXZIaYaFLUiMsdElqhIUuSY2w0CWp\nERa6JDXCQpekRljoktQIC12SGmGhS1IjehZ6kpuTHEny0IplH0/yZJID3eM9mxtTktRLP3vonwUu\nWmP5p6tqR/e4Y7SxJEmD6lnoVfUN4NkxZJEkDWGYY+jXJnmgOyRz6sgSSZI2ZKOFfgPwJmAHcBj4\n5HoDk+xKsphkcWlpaYObkyT1sqFCr6qnq+qlqnoZ+Ayw8zhj91TVQlUtzM3NbTSnJKmHDRV6ku0r\n3r4feGi9sZKk8djSa0CSzwMXAFuTPAH8HnBBkh1AAYeAD29iRklSH3oWelVducbimzYhiyRpCF4p\nKkmNsNAlqREWuiQ1wkKXpEZY6JLUCAtdkhphoUtSIyx0SWqEhS5JjbDQJakRPS/916vP/O7bJx1B\n0ga4hy5JjbDQJakRFrokNcJCl6RGWOiS1AgLXZIaYaFLUiMsdElqhIUuSY2w0CWpERa6JDXCQpek\nRljoktQIC12SGmGhS1IjLHRJaoSFLkmNsNAlqRE9Cz3JzUmOJHloxbLXJ7kryaPd86mbG1OS1Es/\ne+ifBS5atWw3cHdVnQXc3b2XJE1Qz0Kvqm8Az65afBmwt3u9F3jfiHNJkga00WPo26rqcPf6KWDb\niPJIkjZoy7BfoKoqSa23PskuYBfAmWeeOezmXlXmd98+6QiSZshG99CfTrIdoHs+st7AqtpTVQtV\ntTA3N7fBzUmSetlooe8DrupeXwV8ZTRxJEkb1c9pi58H/gl4c5InklwNXA+8K8mjwC907yVJE9Tz\nGHpVXbnOqgtHnEWSNASvFJWkRljoktQIC12SGjH0eeiSZtMkr3M4dP0lE9t2y9xDl6RGWOiS1AgL\nXZIaYaFLUiMsdElqhIUuSY3wtEVpwrxNskbFPXRJaoSFLkmNsNAlqREWuiQ1wkKXpEZY6JLUCAtd\nkhrheeh98DxhSbPAPXRJaoSFLkmNsNAlqREWuiQ1wkKXpEZY6JLUCAtdkhphoUtSIyx0SWqEhS5J\njbDQJakRQ93LJckh4AXgJeBoVS2MIpQkaXCjuDnXO6rqmRF8HUnSEDzkIkmNGLbQC/h6kv1Jdq01\nIMmuJItJFpeWlobcnCRpPcMW+tur6lzgYuCaJOevHlBVe6pqoaoW5ubmhtycJGk9QxV6VT3ZPR8B\nbgN2jiKUJGlwGy70JCclOfnYa+DdwEOjCiZJGswwZ7lsA25LcuzrfK6q7hxJKknSwDZc6FX1OPCW\nEWaRJA3B0xYlqREWuiQ1wkKXpEZY6JLUCAtdkhphoUtSIyx0SWqEhS5JjbDQJakRFrokNcJCl6RG\nWOiS1AgLXZIaYaFLUiMsdElqhIUuSY2w0CWpEcP8CTpJ2pD53bdPZLuHrr9kItsdF/fQJakRFrok\nNcJCl6RGWOiS1AgLXZIaYaFLUiMsdElqxMychz6p81YltWOSPTKOc+DdQ5ekRljoktQIC12SGjFU\noSe5KMm3kzyWZPeoQkmSBrfhQk9yAvBnwMXA2cCVSc4eVTBJ0mCG2UPfCTxWVY9X1Q+AvwEuG00s\nSdKghin004Dvrnj/RLdMkjQBm34eepJdwK7u7YtJvr3Z21zHVuCZCW17I2YtL8xeZvNurlnLC5uY\nOZ8Y6uNv6GfQMIX+JHDGivend8v+j6raA+wZYjsjkWSxqhYmnaNfs5YXZi+zeTfXrOWF2cy80jCH\nXP4FOCvJG5P8EHAFsG80sSRJg9rwHnpVHU1yLfD3wAnAzVX18MiSSZIGMtQx9Kq6A7hjRFk228QP\n+wxo1vLC7GU27+aatbwwm5lfkaqadAZJ0gh46b8kNaKpQk/yI0n+Ocm/JXk4ye+vMeaDSZaSHOge\nvzaJrKsynZDkX5N8dY11P5zkC93tFe5LMj/+hP8v0/HyTuP8HkryYJdncY31SfIn3Rw/kOTcSeRc\nkadX3guSPLdijn93EjlX5Dklya1JvpXkYJK3rVo/bfPbK+9Uze8gZuZ+6H36PvDOqnoxyYnAvUm+\nVlXfXDXuC1V17QTyrec64CDw42usuxr496r6qSRXAJ8Afnmc4dZwvLwwffML8I6qWu/84ouBs7rH\nzwM3dM+TdLy8AP9YVe8dW5rj+2Pgzqq6vDvj7cdWrZ+2+e2VF6ZrfvvW1B56LXuxe3ti95jqXxIk\nOR24BLhxnSGXAXu717cCFybJOLKtpY+8s+gy4C+7fz/fBE5Jsn3SoWZBktcB5wM3AVTVD6rqP1YN\nm5r57TPvzGqq0OGVwwEHgCPAXVV13xrDfrH70e/WJGessX6c/gj4beDldda/couFqjoKPAf8xHii\nralXXpiu+YXlb+pfT7K/u3J5tWm7jUWvvABv6w4tfi3JT48z3CpvBJaAv+gOw92Y5KRVY6ZpfvvJ\nC9MzvwNprtCr6qWq2sHylas7k/zMqiF/B8xX1c8Cd/G/e79jl+S9wJGq2j+pDIPoM+/UzO8Kb6+q\nc1n+0f+aJOdPOlAPvfLeD7yhqt4C/Cnwt+MOuMIW4Fzghqo6B/hPYJpvpd1P3mma34E0V+jHdD9G\n3QNctGr596rq+93bG4GfG3e2Fc4DLk1yiOW7Vb4zyV+vGvPKLRaSbAFeB3xvnCFX6Jl3yuYXgKp6\nsns+AtzG8p1CV+rrNhbj0itvVT1/7NBidy3IiUm2jj3osieAJ1b8JHwry4W50jTNb8+8Uza/A2mq\n0JPMJTmle/2jwLuAb60as/LY3aUs/3JvIqrqd6rq9KqaZ/nWCf9QVb+yatg+4Kru9eXdmIn8XqCf\nvNM0vwBJTkpy8rHXwLuBh1YN2wf8anc2xluB56rq8JijAv3lTfKTx36PkmQny/+PJ/JNvqqeAr6b\n5M3doguBR1YNm5r57SfvNM3voFo7y2U7sDfLf3zjNcAXq+qrSf4AWKyqfcBvJbkUOAo8C3xwYmnX\nsSrvTcBfJXmM5bxXTDTcGqZ8frcBt3X/P7cAn6uqO5P8OkBV/TnLVzu/B3gM+C/gQxPKCv3lvRz4\njSRHgf8GrpjUN/nObwK3dGeMPA58aIrnF3rnnbb57ZtXikpSI5o65CJJr2YWuiQ1wkKXpEZY6JLU\nCAtdkhphoUtSIyx0SWqEhS5JjfgftMywNaaaOAkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10694f240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.952\n",
      "start move from cell 3\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD8CAYAAACRkhiPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEUVJREFUeJzt3X+MZWV9x/H3p6y2lVJBd4oILEtaQoJGkE5WqZagKOVX\noD9o3U1bxWpWDVRtmpi1TbCx/9C01lYx0C1QtFWkUrG0LAixTdDEX7O44PqDinSVXZFdQUGU1q5+\n+8ectcPlzs5wz52Zu33er+Rmznme557nu082+9lz7rlnUlVIktr0EytdgCRp5RgCktQwQ0CSGmYI\nSFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIatWukChlm9enWtXbt2pcuQpAPG1q1bv1VVU0/2fRMZ\nAmvXrmVmZmaly5CkA0aSr43yPi8HSVLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIYZApLUMENAkhpm\nCEhSwybyG8OSHm/tpptWbO4dl56zYnNr6XkmIEkNW/BMIMnVwLnA7qp6btd2HXB8N+RQ4DtVddKQ\n9+4Avgv8ENhbVdNjqluSNAaLuRx0DXAZ8L59DVX1in3bSd4BPLyf97+kqr41aoGSpKWzYAhU1e1J\n1g7rSxLgt4CXjrcsSdJy6PuZwC8DD1TVV+bpL+DWJFuTbNzfgZJsTDKTZGbPnj09y5IkLUbfENgA\nXLuf/hdX1cnAWcBFSU6db2BVba6q6aqanpp60r8XQZI0gpFDIMkq4NeB6+YbU1W7up+7gRuAdaPO\nJ0kavz5nAi8DvlxVO4d1Jjk4ySH7toEzgO095pMkjdmCIZDkWuCTwPFJdiZ5Tde1noFLQUmenWRL\nt3s48IkkdwKfAW6qqlvGV7okqa/F3B20YZ72C4e0fQM4u9u+FzixZ32SpCXkYyM0shYfZbCSf2Zp\nKfjYCElqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1\nzBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhC4ZAkquT7E6yfU7bnyTZlWRb9zp7nvee\nmeTuJPck2TTOwiVJ/S3mTOAa4Mwh7e+sqpO615bBziQHAe8BzgJOADYkOaFPsZKk8VowBKrqduCh\nEY69Drinqu6tqh8AHwTOH+E4kqQl0uczgYuT3NVdLjpsSP+RwH1z9nd2bUMl2ZhkJsnMnj17epQl\nSVqsUUPgcuDngZOA+4F39C2kqjZX1XRVTU9NTfU9nCRpEUYKgap6oKp+WFU/Av6W2Us/g3YBR8/Z\nP6prkyRNiJFCIMkRc3Z/Ddg+ZNhngeOSHJvkqcB64MZR5pMkLY1VCw1Ici1wGrA6yU7gbcBpSU4C\nCtgBvK4b+2zgyqo6u6r2JrkY+ChwEHB1VX1hSf4UkqSRLBgCVbVhSPNV84z9BnD2nP0twBNuH5Uk\nTQa/MSxJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCk\nhhkCktQwQ0CSGmYISFLDDAFJapghIEkNW/A3i2nyrd1000qXIOkA5ZmAJDVswRBIcnWS3Um2z2n7\n8yRfTnJXkhuSHDrPe3ck+XySbUlmxlm4JKm/xZwJXAOcOdB2G/Dcqnoe8B/AW/fz/pdU1UlVNT1a\niZKkpbJgCFTV7cBDA223VtXebvdTwFFLUJskaYmN4zOB3wNunqevgFuTbE2ycQxzSZLGqNfdQUn+\nGNgLvH+eIS+uql1Jfg64LcmXuzOLYcfaCGwEWLNmTZ+yJEmLNPKZQJILgXOB366qGjamqnZ1P3cD\nNwDr5jteVW2uqumqmp6amhq1LEnSkzBSCCQ5E3gLcF5VfX+eMQcnOWTfNnAGsH3YWEnSyljMLaLX\nAp8Ejk+yM8lrgMuAQ5i9xLMtyRXd2Gcn2dK99XDgE0nuBD4D3FRVtyzJn0KSNJIFPxOoqg1Dmq+a\nZ+w3gLO77XuBE3tVJ0laUn5jWJIa5rODdEDyeUnSeHgmIEkNMwQkqWGGgCQ1zBCQpIYZApLUMENA\nkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSp\nYYsKgSRXJ9mdZPuctmckuS3JV7qfh83z3ld1Y76S5FXjKlyS1N9izwSuAc4caNsEfKyqjgM+1u0/\nTpJnAG8DXgCsA942X1hIkpbfokKgqm4HHhpoPh94b7f9XuBXh7z1V4Dbquqhqvo2cBtPDBNJ0grp\n85nA4VV1f7f9TeDwIWOOBO6bs7+za5MkTYCxfDBcVQVUn2Mk2ZhkJsnMnj17xlGWJGkBfULggSRH\nAHQ/dw8Zsws4es7+UV3bE1TV5qqarqrpqampHmVJkharTwjcCOy72+dVwD8PGfNR4Iwkh3UfCJ/R\ntUmSJsBibxG9FvgkcHySnUleA1wKvDzJV4CXdfskmU5yJUBVPQT8KfDZ7vX2rk2SNAFWLWZQVW2Y\np+v0IWNngNfO2b8auHqk6iRJS8pvDEtSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIa\nZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWEj\nh0CS45Nsm/N6JMmbB8acluThOWMu6V+yJGlcVo36xqq6GzgJIMlBwC7ghiFDP15V5446jyRp6Yzr\nctDpwFer6mtjOp4kaRmMKwTWA9fO03dKkjuT3JzkOWOaT5I0Br1DIMlTgfOADw3pvgM4pqpOBN4N\nfGQ/x9mYZCbJzJ49e/qWJUlahHGcCZwF3FFVDwx2VNUjVfVot70FeEqS1cMOUlWbq2q6qqanpqbG\nUJYkaSHjCIENzHMpKMmzkqTbXtfN9+AY5pQkjcHIdwcBJDkYeDnwujltrweoqiuAC4A3JNkLPAas\nr6rqM6ckaXx6hUBVfQ945kDbFXO2LwMu6zOHJGnp9AoBSf//rd1000qXsOx2XHrOSpewbHxshCQ1\nzBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMM\nAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktSw3iGQZEeSzyfZlmRmSH+SvCvJPUnuSnJy3zklSeMx\nrl80/5Kq+tY8fWcBx3WvFwCXdz8lSStsOS4HnQ+8r2Z9Cjg0yRHLMK8kaQHjOBMo4NYkBfxNVW0e\n6D8SuG/O/s6u7f65g5JsBDYCrFmzZgxlLa+1m25a6RIk6Ukbx5nAi6vqZGYv+1yU5NRRDlJVm6tq\nuqqmp6amxlCWJGkhvUOgqnZ1P3cDNwDrBobsAo6es39U1yZJWmG9QiDJwUkO2bcNnAFsHxh2I/DK\n7i6hFwIPV9X9SJJWXN/PBA4Hbkiy71gfqKpbkrweoKquALYAZwP3AN8HXt1zTknSmPQKgaq6Fzhx\nSPsVc7YLuKjPPJKkpeE3hiWpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGG\ngCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGjRwCSY5O8u9Jvpjk\nC0neNGTMaUkeTrKte13Sr1xJ0jj1+UXze4E/rKo7khwCbE1yW1V9cWDcx6vq3B7zSJKWyMhnAlV1\nf1Xd0W1/F/gScOS4CpMkLb2xfCaQZC3wfODTQ7pPSXJnkpuTPGcc80mSxqPP5SAAkvwM8E/Am6vq\nkYHuO4BjqurRJGcDHwGOm+c4G4GNAGvWrOlbliRpEXqdCSR5CrMB8P6q+vBgf1U9UlWPdttbgKck\nWT3sWFW1uaqmq2p6amqqT1mSpEXqc3dQgKuAL1XVX84z5lndOJKs6+Z7cNQ5JUnj1edy0IuA3wU+\nn2Rb1/ZHwBqAqroCuAB4Q5K9wGPA+qqqHnNKksZo5BCoqk8AWWDMZcBlo84hSVpafmNYkhrW++6g\nSbN2000rXYKkA9xK/Tuy49Jzln1OzwQkqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJ\nDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDWsVwgkOTPJ3Unu\nSbJpSP9PJrmu6/90krV95pMkjdfIIZDkIOA9wFnACcCGJCcMDHsN8O2q+gXgncCfjTqfJGn8+pwJ\nrAPuqap7q+oHwAeB8wfGnA+8t9u+Hjg9SXrMKUkaoz4hcCRw35z9nV3b0DFVtRd4GHhmjzklSWO0\naqUL2CfJRmBjt/tokruXaerVwLeWaa5xOJDqPZBqBetdagdSvStSa0a/YL4aOGaUN/YJgV3A0XP2\nj+raho3ZmWQV8HTgwWEHq6rNwOYe9YwkyUxVTS/3vKM6kOo9kGoF611qB1K9B1Kt8ON6147y3j6X\ngz4LHJfk2CRPBdYDNw6MuRF4Vbd9AfBvVVU95pQkjdHIZwJVtTfJxcBHgYOAq6vqC0neDsxU1Y3A\nVcDfJ7kHeIjZoJAkTYhenwlU1RZgy0DbJXO2/wv4zT5zLINlvwTV04FU74FUK1jvUjuQ6j2QaoUe\n9carM5LULh8bIUkNayoEkhyU5HNJ/nVI34VJ9iTZ1r1euxI1zqlnR5LPd7XMDOlPknd1j+S4K8nJ\nK1FnV8tCtZ6W5OE5a3vJsOMslySHJrk+yZeTfCnJKQP9E7O2XT0L1Tsx65vk+Dl1bEvySJI3D4yZ\niPVdZK0Ts7ZdPX+Q5AtJtie5NslPDfQ/6Uf1TMz3BJbJm4AvAT87T/91VXXxMtazkJdU1Xz3Kp8F\nHNe9XgBc3v1cKfurFeDjVXXuslWzf38N3FJVF3R3tj1toH/S1nahemFC1req7gZOgh8/WmYXcMPA\nsIlY30XWChOytkmOBN4InFBVjyX5R2ZvtrlmzrAfP6onyXpmH9Xziv0dt5kzgSRHAecAV650LWNy\nPvC+mvUp4NAkR6x0UZMuydOBU5m9c42q+kFVfWdg2MSs7SLrnVSnA1+tqq8NtE/M+s4xX62TZhXw\n0933rp4GfGOg/0k/qqeZEAD+CngL8KP9jPmN7vT0+iRH72fccijg1iRbu29TD1rMYzuWy0K1ApyS\n5M4kNyd5znIWN+BYYA/wd92lwSuTHDwwZpLWdjH1wuSs71zrgWuHtE/S+u4zX60wIWtbVbuAvwC+\nDtwPPFxVtw4Me9KP6mkiBJKcC+yuqq37GfYvwNqqeh5wG/+XpivlxVV1MrOnzhclOXWF69mfhWq9\nAzimqk4E3g18ZLkLnGMVcDJweVU9H/ge8ITHoE+QxdQ7SesLQHfZ6jzgQytdy0IWqHVi1jbJYcz+\nT/9Y4NnAwUl+p+9xmwgB4EXAeUl2MPu005cm+Ye5A6rqwar67273SuAXl7fEx+tSn6razex1ynUD\nQxbz2I5lsVCtVfVIVT3abW8BnpJk9bIXOmsnsLOqPt3tX8/sP7JzTczasoh6J2x99zkLuKOqHhjS\nN0nrC/updcLW9mXAf1bVnqr6H+DDwC8NjPnx2maBR/Xs00QIVNVbq+qo7tka65l9fMXjEnTgmuR5\nzH6AvCKSHJzkkH3bwBnA9oFhNwKv7O60eCGzp4b3L3Opi6o1ybP2XZdMso7Zv3f7/Yu5VKrqm8B9\nSY7vmk4HvjgwbCLWFhZX7ySt7xwbmP/yysSsb2feWidsbb8OvDDJ07qaTueJ/0496Uf1tHZ30OPk\n8Y+4eGOS84C9zD7i4sIVLO1w4Ibu794q4ANVdUuS1wNU1RXMflP7bOAe4PvAqye41guANyTZCzwG\nrF/hZ0j9PvD+7jLAvcCrJ3Rt91mo3ola3+4/Ay8HXjenbSLXdxG1TszaVtWnk1zP7CWqvcDngM3p\n+agevzEsSQ1r4nKQJGk4Q0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIb9L0G7bjTXTCIb\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106626358>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.1\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "for start in [1,2,3]:\n",
    "    Run_total = np.zeros(100)\n",
    "    N_move    = 50\n",
    "    print(\"start move from cell\", start)\n",
    "    for each_run in range(100):\n",
    "        move = np.zeros(N_move)\n",
    "        for each_move in range(N_move):\n",
    "            move[each_move] = steps(start)\n",
    "        Run_total[each_run] = np.mean(move)\n",
    "    plt.hist(Run_total,bins = 10)\n",
    "    plt.show()\n",
    "    print(Run_total.mean())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}