notebooks/book2/12/random_walk_integers.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a2d4b8ea",
   "metadata": {},
   "source": [
    "# Symmetric random walk on uniform distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "85542952",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example from section 29.4 & 29.6 (Fig 29.14 & 29.15) of https://www.inference.org.uk/itprnn/book.pdf\n",
    "try:\n",
    "    import probml_utils as pml\n",
    "except ModuleNotFoundError:\n",
    "    %pip install -qq git+https://github.com/probml/probml-utils.git\n",
    "    import probml_utils as pml\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import jax.numpy as jnp\n",
    "import seaborn as sns\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9d56e221",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/patel_zeel/miniconda3/envs/py3713_2/lib/python3.7/site-packages/probml_utils/plotting.py:26: UserWarning: LATEXIFY environment variable not set, not latexifying\n",
      "  warnings.warn(\"LATEXIFY environment variable not set, not latexifying\")\n"
     ]
    }
   ],
   "source": [
    "pml.latexify(width_scale_factor=2.2, fig_height=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce4f9876",
   "metadata": {},
   "source": [
    "We genereate samples from following distribution\n",
    "$$\n",
    "P(x)= \\begin{cases}1 / 21 & x \\in\\{0,1,2, \\ldots, 20\\} \\\\ 0 & \\text { otherwise. }\\end{cases}\n",
    "$$\n",
    "The proposal distribution is\n",
    "$$\n",
    "Q\\left(x^{\\prime} ; x\\right)= \\begin{cases}1 / 2 & x^{\\prime}=x \\pm 1, x \\in\\{1,2, \\ldots, 19\\} \\\\ 1 / 2 & x' = \\{0,1\\}, x = 0  \\\\1 / 2 &x' = \\{20, 19\\}, x=20 \\\\ 0 & \\text { otherwise. }\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "60a0a67c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create transition matrix\n",
    "prob_vect = jnp.repeat(0.5, 20)\n",
    "transition_mat = jnp.diag(prob_vect, -1) + jnp.diag(prob_vect, 1)\n",
    "transition_mat = transition_mat.at[0, 0].set(0.5)\n",
    "transition_mat = transition_mat.at[20, 20].set(0.5)\n",
    "\n",
    "states = [0, 1, 2, 3, 10, 100, 200, 400]\n",
    "MARKER_SIZE = 2 if pml.is_latexify_enabled() else 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "3e2e1364",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_q_at_state_no(state_no, p_initial, transition_mat):\n",
    "    transition_mat_at_state_no = jnp.linalg.matrix_power(transition_mat, state_no)\n",
    "    q = jnp.matmul(transition_mat_at_state_no, p_initial)  # q = p_initial * T^n = transpose(T)^n * p_initial\n",
    "    return q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5550e614",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_stem(ax):\n",
    "    markerline, stemlines, baseline = ax.stem(q, markerfmt=\"ro\", linefmt=\"blue\")\n",
    "    ax.set_yticks([])\n",
    "    ax.set_ylim(0, max(q) * 1.13)\n",
    "    ax.set_ylabel(\"$p_{\" + str(state_no) + \"}(x)$\", fontdict={\"rotation\": 0}, labelpad=13)\n",
    "    baseline.set_linewidth(0)\n",
    "    stemlines.set_linewidth(1)\n",
    "    markerline.set_markersize(MARKER_SIZE)\n",
    "    sns.despine()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8883cc66",
   "metadata": {},
   "source": [
    "## $x_0 = 10$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "74dc9505",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DeviceArray([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
       "             0], dtype=int32)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_0 = 10  # initial state\n",
    "p_initial = [0] * 21\n",
    "p_initial[x_0] = 1\n",
    "p_initial = jnp.array(p_initial)\n",
    "p_initial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "dc0904ea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEYCAYAAAC0tfaFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAiuUlEQVR4nO3df7QcZZ3n8fcnNwmBiE4E9GRAVi/jEC9kA1zgyMp44BxwL2LAYWYHhgw64pFl/MEgzHpgdDeL4Lo6B5kDsjjLwgqKmB0mBhHMoCwQBQ4/VxwSfg3ZZMEfyEVwQMWQ5Lt/1HNNc+m+t7u6bz+3uj6vc/p091P19Pep6upvVz9dVY8iAjMzq5c5uRtgZmb95+RvZlZDTv5mZjXk5G9mVkNO/mZmNeTkb2ZWQ07+ZmY15ORvM0bSeklHlJ3eMN8mSUf1rmXlSfqypAvS42nb3+4ymvWbk791pJNEHBH7RcRtreo1Tu9Bu06WdJ+kFyX9RNK3JR3ei9duZXL7Z3oZc5H00bRufyPpy02mv17SNyT9UtJmSSdnaKZ1aG7uBph1S9JZwDnA6cA/AluAMeB44PsZmzYofgxcAPxbYOcm0y+lWOdvBA4AbpT0YESs71sLrWPe87fS0p7uX0n6oaRfSFolacGk6UdJ+gqwN3BD2jP/ROP09PgcSU9IekHSBkl/2GYbXgd8GvhIRKyOiF9GxMsRcUNE/Ic0z9sk3Sbp+dQNc1wHy3CgpAdSu1YBr1q+9LidZSzdjibL/XlJaxqe/42kWyTNb2e9dSKt1zXAs03asRD4I+A/RsSLEfF94JvAKb1uh/WWk791608o9rLfAvxr4M8nzxARpwD/D1geEa+JiM83eZ0ngD8AXgecB3xV0uI24h9GkZC/0WyipHnADcDNwBuAjwHXSNp3umVIiXQN8BXg9cDfUyS6V5luGbtpRwufA45MX06np3onRMSWKepMtOVb6Quo2e1b09Wf5PeBrRHxWEPZg8B+Hb6O9ZmTv3Xr4oj4cUT8nCK5HVDmRSLi79PrbI+IVcDjwKFtVN0NGI+IrS2mvx14DfBfI2JLRPxv4FvAn7axDG8H5gF/m35NXAfc2/HCdd+OV4mIZ4GLgKuAc4F3R8QvACR9TtL3JH0lfelMrvueiPidFrf3dLhcrwH+ZVLZL4BdO3wd6zMnf+vWTxse/4oiGXRM0vsk/WBiDxTYH9i9jarPArtLavX/1e8CT0bE9oayzcCeDc9bLcPvAj+KV176dnMbbep1O1r5P8BS4NyIeBJA0jJgz4j4A+AR4I9LtrddLwKvnVT2WuCFGY5rXXLyt35pee1wSf8KuBz4KLBbRPwO8BCgNl73LuA3wHtbTP8x8CZJjdv63sCP2njtnwB7Smpsx95TzD/V9dG7acerSFoKXEax539qw6R/Q9G1BLAWeEeTut9O/0s0u327w6Y8BsyV9NaGsmWA/+yd5Zz8rV+eBoZbTFtIkTifAZD0AYo9/2mlro7/BFwq6b2SdpE0T9Ixkj4P3E2xF/2JVH4EsBz4ehsvfxewFTgj1T2BqbuiplrGbtrxCpL2pOgWOh34MLBUO84lWMSObphfUPxX8QoRcUz6X6LZ7Zgm8eamP5+HgCFJCyZ+aUXEL4HVwKclLZT0DoqjrL7S6XJZfzn5W798FvhU6tb5q8YJEbEBuJAi2T5N0ZVxR7svHBEXAmcBn6L4AnmS4lfEmvQH6HLgGGAc+G/A+yLikTZedwtwAsUfrz8HTqRIdGWWsXQ7Gkl6LXAT8IWI+GZE/Ar4G+AzaZbn2dEN87rU7m59Cvg1xeG0f5Yef6ph+ocpDgH9GXAt8Bc+zHP2k0fyMhsckg4AzoqI90n6a+D/RsS1mZtls5D3/M0GSET8AHha0vcoDrf8h7wtstnKe/5mZjXkPX8zsxpy8jczqyEnfzOzGqrEVT3HxsZi7dq1uZthZlY1LU+U7Muev6SdJd0uaajF9PmS1rU6RX98fHxmG2hWxsaNsN9+MHducb9xY+4WmbWtX90+pwKrI2Jbs4npBJhbKE6iMauG5cvhkUdg27bifvny3C0ya1vp5C/p2nTN8XtUjN5z7BSzrwCub6h7q6Sj0+MLJF1CcencFWXbY9Z3jz4K29N12rZvL56bVUQ3e/7LgI0RcShF0l7ZbKZ0TfThiNjUULwS+KSkFcCBwJkUF/I6pIv2mPXXvvvCnPQRmjOneG5WEaWSf7rI0x4Ug24AbAAWpQs7XSXp8pTYobgs7/ON9SNiHcUfEWcBJ0XEttQltEWSrwNu1XDDDbBkCVsZgiVLiudmFVH2aJ/9gccj4qX0/CCK0XtOAK6LiBvSkHfXUFwE6hXD0aXL0S4Gno2Ixut+7wS8hFkVDA/D+vXME/gyZlY1Zbt9lgF7p0u7LqT4BXARsBfFFRUBtgFExHOky8ACpKH5rqG47OuLksZS+cSITC+XXRgzM2tPN8l/NcU1yu8FLouIO4CnKL4AJr/2zcDhknZJ9c6OiIeB89nxX8GRwI0l22NmZh0o2+2zDDgtIs6YVL4a+GI68qexA/RS4OMR8V2KAbeB3/b9Tzw/meJ64WZmNsPKJv99KAbYfoU0qs8HmpQ/kA7vHGp2rH86ImhNRDxWsj1mZtaBUsk/Ivaafq5X1blyimlbgKvLtMXMzDrnC7uZmdWQk7+ZWQ05+ZuZ1ZCTv5lZDTn5m5nVkJO/mVkNOfmbmdWQk7+ZWQ05+ZuZ1dDgJv9ux1ftpn7d6la13TnH4K3iMnsbqc4ytyMiZv1tdHQ0OjYyEjFnTgQU9yMj/atft7pVbXe3y5xAiUpVXGZvI9VZ5h1a5lVFRG+/TWaApGeAzZ3UGYXRyWX3w/39qF+3ujlj51zmblRxmb2NdFY3d+xkPCLGmk2oRPI3M7PeGtw+fzMza8nJ38yshpz8zcxqqOxIXn01NjYWa9euzd0MM7OqUasJfdnzl7SzpNslDbWYPl/SOklNv4zGx8dntoFWT1U9zt+sB/rV7XMqsDqajN8Lvx3G8RbgxD61xwyWL4dHHoFt24r75cvrEduMLpK/pGslrZJ0j6TNko6dYvYVwPUNdW+VdHR6fIGkS4A1aT6z/nj0Udi+vXi8fXvxvA6xzehuz38ZsDEiDqVI2iubzSRpPjAcEZsailcCn5S0AjgQOBN4CDiki/aYdWbffWFO+gjMmVM8r0NsM0omf0kLgD2A81LRBmCRpGFJV0i6rmH23YHnG+tHxDqKPyLOAk6KiG2pS2iLpF3LtMmsYzfcAEuWsJUhWLKkeF6H2GaUP9pnf+DxiHgpPT8IeDAiNgIfnJT8fw0saKwsaSmwGHg2Il5omLQT8BJm/TA8DOvXM08Q62sU24zy3T7LgL0lLZC0kOIXwEXNZoyI54Ch9GsBSYuBa4DjgRcljaXy3SiuQ/FyyTaZmVmbukn+q4G7gXuByyLijinmvxk4XNIuqd7ZEfEwcD47/is4ErixZHvMzKwD3ST/SyNiWUSMRMRXodh7l/Ql4EBJ5zbMfynw/oj4VUQcFhHfgaLvPyIOS/OcDPxdyfaYmVkHyvb57wM8PrkwIp4FTm9S/kA6vHOo2bH+6YigNRHxWMn2mJlZB0ol/4jYq0SdK6eYtgW4ukxbzMysc76wm5lZDTn5m5nVkJO/mVkNOfmbmdWQk7+ZWQ05+ZuZ1ZCTv5lZDTn5m5nVkJO/mVkNDW7y73aM1G7q161uldvdjVyxq7quq9juqi5zOyJi1t9GR0ejYyMjEXPmREBxPzLSv/p1q1vldkdRtZRcsau6rqvY7qou8w4t86oiorffJjNA0jPA5k7qjMLo5LL74f5+1K9b3Zyxu213N3LFruq6rmK7q7rMDcYjYqzZhEokfzMz663B7fM3M7OWnPzNzGrIyd/MrIb6kvwl7SzpdklDLabPl7ROUtPBZcbGxgLwzTfffPOts1tL/drzPxVYHU2GcAQmRvK6BTix2fTx8fEZbJpZSTnPMTDrUunkL+laSask3SNps6Rjp5h9BXB9Q91bJR2dHl8g6RJgTZrPrBqWL4dHHoFt24r75ctzt8isbWUHcAdYBlwfESdKOhz4AnDj5JnS4OzDEbGpoXgl8GlJbwAOBI5L5Yd00R6z/nr0Udi+vXi8fXvx3KwiSu35S1oA7AGcl4o2AIskvVfS5ekXwbvStN2B5xvrR8Q6QMBZwEkRsS11CW2RtGuZNpn13b77wpz0EZozp3huVhFlu332Bx6PiJfS84OAByNiTUR8CDidHf33vwYWNFaWtBRYDGyJiBcaJu0EvIRZFdxwAyxZwlaGYMmS4rlZRZRN/suAvSUtkLSQ4hfARQ3TPwVcChARzwFD6dcCkhYD1wDHAy9KGkvlu1GcivxyyTaZ9dfwMKxfzzy2wvr1xXOziugm+a8G7gbuBS6LiDtU+Bzw7Yh4oGH+m4HDJe2S6p0dEQ8D51P0/wMcSZP/DMzMrPfK/uG7DDgtIs6YVP4x4CjgdZJ+LyK+lMovBT4eEd8FDpuYOfX9Tzw/GTinZHvMzKwDZZP/PsDjkwsj4mLg4iblD6TDO4eaHeufjghaExGPlWyPmZl1oFTyj4i9StS5coppW4Cry7TFzMw652v7WLXV8SzbOi6z9ZyTv1VbHc+yreMyW885+Vu11fEs2zous/Wck79VWx3Psq3jMlvPOflbtdXxLNs6LrP1XCXG8D344IPjvvvuy90Mm8UkyLUp54qdc5mtMtRqgvf8zcxqyMnfzKyGnPzNzGrIyd/MrIYGN/l3exZkN/XrVjd37LJynilbxWWu6jZSx2VuR0TM+tvo6Gh0bGQkYs6cCCjuR0b6V79udXPHjqJqx3oQN2fsvset6jZSx2XeoWVercShnpKeATZ3UmcURieX3Q/396N+3ermjl1Wrrg5Y9dxG6njMjcYj4ixZhMqkfzNzKy3BrfP38zMWnLyNzOrISd/M7Ma6kvyl7SzpNslDbWYPl/SOklNRxYbGxsLwDffenvbuDHYb79g7tzifuPGesT2rU63lvq1538qsDqajN8LTAzjeAtwYrPp4+PjM9g0q62cg6J4QBbLrHTyl3StpFWS7pG0WdKxU8y+Ari+oe6tko5Ojy+QdAmwJs1n1h85B0XxgCyWWTd7/suAjRFxKEXSXtlsJknzgeGI2NRQvBL4pKQVwIHAmcBDwCFdtMesMzkHRfGALJZZqeQvaQGwB3BeKtoALJL0NklfknSdpL9I03YHnm+sHxHrKK4zfRZwUkRsS11CWyTtWqZNZh3LOSiKB2SxzJr+wdqG/YHHI+Kl9Pwg4MGIeBg4XdIc4GrgMuDXwILGypKWAouBZyPihYZJOwEvYdYPw8Owfj3zBLG+RrHNKN/tswzYW9ICSQspfgFcBCDpOOBG4CaAiHgOGEq/FpC0GLgGOB54UdJYKt+N4lTkl7tYHjMza0M3yX81cDdwL3BZRNwBEBHfjIhjeOWftzcDh0vaJdU7O/1KOJ8d/xUcSfGlYWZmM6xst88y4LSIOKOxUNIRwAkU3Tc3NUy6FPh4RHwXOGyiMPX9Tzw/GTinZHvMzKwDZZP/PsDjkwsj4jbgtiblD6TDO4eaHeufjghaExGPlWyPmZl1oFTyj4i9StS5coppWyj+IDYzsz7wtX3MzGrIyd/MrIac/M3MasjJ3/LLOZB63XhdW+Lkb/n5Cpf943VtiZO/5ecrXPaP17UlTv6Wn69w2T9e15Y4+Vt+vsJl/3hdW6KIKUf6mhUOPvjguO+++3I3w2aYBDk2x1xxc8bOuczWV2o1wXv+ZmY1NLjJv9tD2rqpX7e6vahfVs5DF+u2zP5MVWeZ2xERs/42OjoaHRsZiZgzJwKK+5GR/tWvW91e1I+iasdyxc0Zu6pxq7htV3WZd2iZVyvR5y/pGWBzJ3VGYXRy2f1wfz/q161uL+qXlStuzthVjVvFbbuqy9xgPCLGmk2oRPI3M7PeGtw+fzMza8nJ38yshsqO5NVXY2NjsXbt2tzNMDOrmvzH+UvaWdLtkoammGe+pHWSXvGlND4+PvMNNOuUr5BpFdbPbp9TgdXRZAzfCVEM53gLcGLfWmVWlq+QaRXWUfKXdK2kVZLukbRZ0rEdVF8BXN/wWrdKOjo9vkDSJWnSmjSv2ezmK2RahXW6578M2BgRh1Ik6JXtVJI0HxiOiE0NxSuBT0paARwInJnKHwIO6bBdVlc5u15yXSHT3U3WA20nf0kLgD2A81LRBmCRpGFJV0i6Ls23UNJVki5PiR1gd+D5xteLiHUUf0acBZw00R2U7rdI2rWL5bK6yNn1kusKme5ush7oZM9/f+DxiHgpPT8IeDAiNkbEBxvmOwG4LiI+BByXyn4NLGh8MUlLgcXAloh4YVKsnYCXMJtOzq6X4WFYv555bIX164vn/eDuJuuBTpL/MmBvSQskLaT4BXBRk/n2Ap5Mjyf25p8DhtKvByQtBq4BjgdelPTb048l7UZxSvLLnS6M1VAdByep4zJbz3Wa/FcDdwP3ApdFxB1N5nuK4gtg8uvfDBwuaZf0OmdHxMPA+bzyv4MjgRs7aJfVWR0HJ6njMlvPtX1tH0m3A6dFxKOTyncDPgMcDfwP4GLgixTdNt+PiGvSfAcBH4+IU6aJsxo4JyIemyjzYC42HQ/IYtZUy5O8OjnDdx/g8cmFEfEscPqk4g80me+BdHjnUKtj/dNRQWsaE7+ZmfVe28k/Ivaafq5pX+PKaaZvAa7uNo6ZmU3NF3az3vCx54PP7/FAcfK33vCx54PP7/FAcfK33vCx54PP7/FAcfK33vCx54PP7/FAcfK33vCx54PP7/FAqcQYvj7OvzrqeMx73ZbZ5xdUSv7BXKwPfDSGDSpv2z3n5D9IfDSGDSpv2z3n5D9IfDSGDSpv2z03uMk/58/EXLF9NIYNqjoOnDPTsSNi1t9GR0ejYyMjEXPmREBxPzLSWf0nnijqDA0V90880Z/Y3cRNdV+mfN1ScRtAiUo9iJ0rbs7YlYrbbexc23ZV88gOLfNq9sTezq1U8h8aKhZv4jY01Fn9blZ8N7F78IaX+nD2ZkPLFtvLPMvj5oyd67OcO3ahZV6txKGekp4BNndSZ3/Yb6eG0cN+Ay89BOvbrT8Ko5PL7of7Zzp2N3G7kStuzthe5v7FzRk712c5d+xkPCLGmk2oRPI3M7PeGtw/fM3MrCUnfzOzGupkJK9sxsbGYu3atbmbYWZWNbPr8g6SdpZ0u6ShFtPnS1onaS7A+Ph4fxto1THIx2HPtrg2UHJ1+5wKrI4WY/lGMZzjLcCJfW2VlZcrIeU87T9X7Fxx/aUzUNo62kfStRRfFG8B3gh8OCJuLB1UuhM4OSI2SboV+C8R8R1JFwCvi4iPSVoGfDYi3u2relbAfvsViWj79uIMzCVLYH2nR6WVMHdukQQnDA3B1q0zHzdn7Fxxc73H1o2uu32WARsj4lBgBbCydEuk+cBwRGxKRSuBT0paARwInJnKHwIOKRunlnLumeW69krOS1rkip0rbs7r6/hXR89Nm/wlLQD2AM5LRRuARZKGJV0h6bo030JJV0m6PCXypmXA7sDzE68fEesovp3OAk6a6ApK91sk7dqTJa2DnF0guRJSzgFGcsXOFTfnF62v6tl7U53+m7qEDgbubHh+FHBdw/Pr0v0pwPL0eNUUZYuATQ31lwKPAXc1iT0OzCt1eYecenTNmI715nTwcrq59koPlL7kQIVj9z1uzvc417ad67PcOy1zezvdPsuAvSUtkLSQ4hfARU3m2wt4Mj3e1qosIp4DhtLrLQauAY4HXpT029OQJe1GcWryy220cXbJtZeSc89seBjWr2ceW4t+4OHh/sW2/sj5Hufatgf4F0e7yX81cDdwL3BZRNzRZL6nKJJ94+s2KwO4GXhXet2zI+Jh4Hxe+V/CkUDpP5WzytU36jFWbVDl2rYHeByBdpP/pRGxLCJGIuKrUOyZS/oScKCkcykS+R9JugyYeGealQFcCvy7iDgsIr4DRd9/RBzWMM/JwN+VXrKcfxDl2kvx3rcNqlzbds5f0zOcw6Y91FPSU8DeEbG9p4GlU4Grosmx/umIoJMi4mooOYB7zsPSNm6E5cvZuuFR5o7sW+yl9DERezDzesSuW9wssXN+lnuTw1oe6lmJq3qWSv7dHgud3nQefbT4ti/xpvvDOfhxc8auW9ycsUvH7SaP9OZ8jtl1eYe+6Pbn2gD/0WNmfdJNHpnhLqfBTf7d/kGU648en8xi1ltVPflxhv/kHtxun6T0z7Ue9LeVip0rbo9U7md5hWPXLW7p2Dk/U/k/zzXs9ulW3Q4t8y8Om2m5trGch2vO4sOvnfxbqduhZTn/4/AXT//kXNc++XFWHX7t5D/b1O0XB9TzEsV1vAS2T36cVdznP4P1K1W3R+dFlIrdg0PaKrfMufqCc61ryN7/natu5tju87dp5Nw7quMliut4CWzvgc8qTv5WyNk3WcdLFNfxEtizuP+7jpz8Lb9cSaGOYwE4AVviPv8ZrF+3ujlje5mrUTdn7DouM1W/to+kZ4DNudthZlYx4xEx1mxCJZK/mZn1lvv8zcxqyMnfzKyGnPzNzGooS/KXtLOk2yUNtZg+X9I6SXMBxsbGAvDNN998862zW0u59vxPBVY3G8IRICK2ALcAJwKMj4/3sWlWSq5r1dTx+jx1i2szIyKmvQHXAquAeygOuTy2nXpTvN6dwJvT41uBo9PjC4BL0uNlwE0RwejoaFgbnngiYmQkYmiouH/iif7FHhmJmDMnAor7kZHBjpszdt3iRuTdtqutdR6eamLsSNYbgM+mx4cD97RTr8VrzQd+2vD8ncBtwArgRmAolQ8Bz0QVk3+uDTXnh3NoqIg7cRsaGuy4OWPXLW5Evm27+l86LXPxtN0+khYAewDnpaINwCJJ75V0uaRVkt4laaGkq1LZilT3VWXA7sDzDb881lGchXYWcFKkrqB0v0XSrtP/fmmijtctz3mhslzXqqnj9XnqFhfybduDPM7FVN8MUeyBHwzc2fD8KOC6hueLgCuAU4DlqWxVum9WtgjY1FB/KfAYcFeT2OPAvFJ7/t3uKXTzjZ9rD2kW/Cx/mT7vIeWKmzN23eJG5Nu2u/0sd5NHerPM5bt9gA8CTwELgIXAHcA7GqZfCBwEnAsckMq+lu5fVZYeP5lebzHwQ+BtwHeAsYZ5dgMeibLdPt2+ad2s+G7qdrOx5PxwJtD3kFnj5oxdq7jdbNs5E3A39XuzE9lV8r843R6k6PL5s1Qu4HPAUbFjL/896fHXW5Wlx1cAxwF3sePP3nc27v0DfwxcGGWTf7dvWjcrvpsNtQff9qU+nD3q28wV28m/TTnXda7YPdgZK71D1U0emQV7/rcD+zYpPwO4H/gScHr6VfA/gcuAFWmeV5Wl8oOAr0wTdzXw+1E2+Xf7puVKwj34tu/7B2QWxPYXXptyreucsXN9piLyfvEUukr+TwFzppuv0xvFsf5DLabNB9438bybo3263UvpZsVX6sPZo/8pKvXhrOMXXs5EWOH3uXJ5pKF6q9u0R/tExF4RsX26+ToVEVfGFCd5RcTVvY7ZkboNMFLHo2Y8jGM93mePXtaUr+0z29TtSydnbCfCerzPszgB51SJ6/l7JK9q1M0Zu1TdjRth+XK2bniUuSP7FsmoRGLIFbtS63oWxK7jMlP1kbyc/KtRN2dsL3M16uaMXcdlZork724fM7MacvI3M6shJ38zsxpy8jczqyEnfzOzGnLyNzOrISd/M7MacvI3M6shJ38zsxqqxBm+kp6hGDjezMzaNx4RY80mVCL5m5lZb7nbx8yshpz8zcxqyMnfzKyGsiR/STtLul3SUIvp8yWtkzQXYGxsLADfprtt3Bjst18wd25xv3Hj4Mf2Mg9+3Nyxq31rbaoxHmfqBnwE+Mtp5llJGvS9mwHcSw+Q3U39XHV7MFh06fWVK7aXuRpxc8bO9XnMHbvQOsdONTF2JOJrgVXAPRSHXB7bTr0pXu9O4M3p8a3A0enxBcAl6fEy4KYom/y73VC7qZ+rbjcDZHe7vnLF9jJXI27O2Lk+j7ljF7pO/huAz6bHhwP3tFOvxWvNB37a8PydwG3ACuBGYCiVDwHPRNnk3+2G2k39XHX9wfYyz9a4OWPn+jzmjl1omYun7fOXtADYAzgvFW0AFqVpCyXdJ+k96fFVki6XtKJh+ivKgN2B5xu6ndZRDDV2FnBSRGxL5duALZJ2na6NTXU7QHY39XPV7WaA7G7XV67YXuZqxM0ZO9fnMXfs6Uz1zRDFHvjBwJ0Nz48CrkuPPw18AngPcAqwPJWvSvfNyhYBmxpebynwGHBXk9jjwLxu+vxfpru+ulL1c9VNoOMqPYmbM7aXeZbHzRk75+cxcy6IKXL7tGf4SvogxV7/71F0xdycEv4uwG7AgpSklwLfjogfSPpaRJws6dzJZek1nwTemr4I/hE4EbgYuDAi1qZ5dgPuiIglHsC9GnVzxvYyV6Nuzth1XGa6HMB9GbAauBu4F7gsIu4AjgDeDpwMfAj4MbDXpNd9qkkZFF8g70qve3ZEPAycT3GEz4QjKf4DMDOzHpvbxjzLgNMi4ozGwoj4JICkP6fY878V+KKkY4GJDrnVTcoALgU+HhGHNbzeOuCwhnlOBs7paGnMzKwt7ST/fYDHW02MiC83PP3ApGm/nFyWyh+QdKukoUh/8DaSNB9YExGPtdE+MzPr0LTJPyL2mm6eMiLiyimmbQGunom4Zmbma/uYmdWSk7+ZWQ05+ZuZ1ZCTv5lZDTn5m5nVkJO/mVkNOfmbmdWQk7+ZWQ05+ZuZ1ZCTv5lZDTn5m5nVkJO/mVkNTTuYy2wg6RmKgePL2J3iktPWHq+vznmddcbrqzPdrK/xiBhrNqESyb8bku6LiINzt6MqvL4653XWGa+vzszU+nK3j5lZDTn5m5nVUB2S/3/P3YCK8frqnNdZZ7y+OjMj62vg+/zNzOzV6rDnb2Zmkzj5m5nV0EAnf0ljkh6V9M+SzsndntlO0iZJ/yTpB5Luy92e2UbSlZJ+JumhhrLXS/qOpMfT/aKcbZxNWqyv/yzpR2kb+4Gkd+ds42wi6U2SbpW0QdJ6SX+ZymdkGxvY5C9pCLgUOAYYAf5U0kjeVlXCkRFxgI/DburLwOQTZs4BbomItwK3pOdW+DKvXl8AF6Vt7ICIuKnPbZrNtgJnR8QI8HbgIylnzcg2NrDJHzgU+OeI2BgRW4CvA8dnbpNVWESsA34+qfh44Kr0+Crgvf1s02zWYn1ZCxHxk4h4ID1+AXgY2JMZ2sYGOfnvCTzZ8PypVGatBXCzpPslnZa7MRXxxoj4SXr8U+CNORtTER+V9MPULeRusiYkvRk4ELibGdrGBjn5W+cOj4iDKLrKPiLpnbkbVCVRHDftY6endhmwD3AA8BPgwqytmYUkvQb4B+DMiPiXxmm93MYGOfn/CHhTw/O9Upm1EBE/Svc/A75B0XVmU3ta0mKAdP+zzO2Z1SLi6YjYFhHbgcvxNvYKkuZRJP5rImJ1Kp6RbWyQk/+9wFslvUXSfOAk4JuZ2zRrSVooadeJx8C7gIemrmUU29T70+P3A9dnbMusN5HEkj/E29hvSRJwBfBwRHyhYdKMbGMDfYZvOozsb4Eh4MqI+EzeFs1ekoYp9vYB5gJf8/p6JUnXAkdQXGL3aWAlsAb4X8DeFJcd/5OI8J+ctFxfR1B0+QSwCfj3Df3ZtSbpcOB7wD8B21PxX1P0+/d8Gxvo5G9mZs0NcrePmZm14ORvZlZDTv5mZjXk5G9mVkNO/mZmNeTkb2ZWQ07+ZmY15ORvVkK67vrR6fEFki7J3SazTszN3QCziloJfFrSGyiuvnhc5vaYdcRn+JqVJOl24DXAEen662aV4W4fsxIkLQUWA1uc+K2KnPzNOpSuTHkNxQhLL0pqNlSh2azm5G/WAUm7AKspxlp9GDifov/frFLc529mVkPe8zczqyEnfzOzGnLyNzOrISd/M7MacvI3M6shJ38zsxpy8jczq6H/D0z67O9ousR8AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(len(states), 1)\n",
    "axs[0].set_title(f\"Initial Condition $x_0 = {x_0}$\")\n",
    "for ind, state_no in enumerate(states):\n",
    "    q = calculate_q_at_state_no(\n",
    "        state_no, p_initial, transition_mat\n",
    "    )  # q is proposed probabilities for p(x) on particular markov state\n",
    "    plot_stem(axs[ind])\n",
    "    if state_no != states[-1]:\n",
    "        axs[ind].set_xticks([])\n",
    "    else:\n",
    "        axs[ind].set_xticks(range(0, 21, 5))\n",
    "        axs[ind].set_xlabel(\"$x$\")\n",
    "\n",
    "pml.savefig(f\"randomWalk{x_0}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "e7ed5134",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.99999994 0.047619047619047616 [0.04665894 0.04683788 0.04656532 0.04708623 0.04656531 0.04738191\n",
      " 0.04665894 0.04769868 0.04683788 0.04800836 0.04708624 0.04828346\n",
      " 0.04738192 0.04849952 0.04769868 0.04863734 0.04800837 0.04868469\n",
      " 0.04828347 0.04863734 0.04849952]\n"
     ]
    }
   ],
   "source": [
    "print(sum(q), 1 / 21, q)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41f32c3c",
   "metadata": {},
   "source": [
    "## $x_0 = 17$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "737d9513",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DeviceArray([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,\n",
       "             0], dtype=int32)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_0 = 17  # initial state\n",
    "p_initial = [0] * 21\n",
    "p_initial[x_0] = 1\n",
    "p_initial = jnp.array(p_initial)\n",
    "p_initial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "d7b692ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEYCAYAAAC0tfaFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAiTUlEQVR4nO3df5RcZZ3n8fenOwmBwDgZQE+GyGozLrEhG6CBIzOMC+eA24gBR2cHhgw64pFlR2UQZj0weDaL4Lo6B5kDsnGWhRUUMTtMJohgBmWBKHD4EVYcEggZMsmCP5BGcEDFkOS7f9ynodKp6u66t7pu3bqf1zn31K3n3qee57m3+1u3nrr1PIoIzMysXgbKroCZmXWfg7+ZWQ05+JuZ1ZCDv5lZDTn4m5nVkIO/mVkNOfibmdWQg7/NGEnrJR2Xd3vDflskndC5muUn6cuSLkvrU9Z/um006zYHf2tLO4E4Ig6JiLtb5Wvc3oF6nSHpYUkvS/qxpG9JOrYTr93KxPrPdBvLIulj6dj+WtKXm2x/ecKyQ9JVJVTV2jCr7AqYFSXpfOBC4BzgH4BtwChwKvC9EqvWL34EXAb8O2DPiRsjYu/xdUl7Az8B/rZrtbNcfOVvuaUr3b+Q9ANJP5e0UtLcCdtPkPQV4EDg1nRl+MnG7Wn9QklPSXpJ0gZJfzDNOrwB+DTw0YhYFRG/iIhXI+LWiPhPaZ+3S7pb0oupG+aUNtpwuKRHUr1WAru1L61Pp42569Gk3Z+XtLrh+V9JulPSnOkct3ak47oaeH4au78f+Cnw3U7XwzrLwd+K+iOyq+y3Av8G+NOJO0TEmcD/A5ZGxN4R8fkmr/MU8PvAG4BLgK9KWjCN8o8hC8h/32yjpNnArcAdwBuBjwM3Sjp4qjakQLoa+ArwW2RXs+9vVs5UbSxSjxY+Bxyf3pzOSfneFxHbJskzXpdvpjegZss3p8o/hQ8CN4QHDet5Dv5W1JUR8aOI+BlZcDssz4tExN+m19kZESuBTcDR08i6LzAWEdtbbH8HsDfw3yJiW0T8H+CbwB9Pow3vAGYDf50+TdwMPNR244rXYzcR8TxwBXA9cBHw7oj4OYCkz0n6rqSvpDediXnfExG/2WJ5T872IelfAf821cl6nIO/FfWThvVfkgW4tkn6gKTvj1+BAocC+00j6/PAfpJafX/128DTEbGzIW0rcEDD81Zt+G3ghxOuYrdOo06drkcr/xdYDFwUEU8DSFoCHBARvw88AfxhzvrmcSbwvYj45y6WaTk5+Fu3tOwGSFeM1wAfA/aNiN8EHgM0jde9H/g18N4W238EvFlS49/6gcAPp/HaPwYOkNRYjwMn2X+yro4i9diNpMXACrKr7LMaNv0uWdcSwBrg95rk/VaTO3TGl2/lqU/yAXzVXxkO/tYtzwJDLbbNIwuczwFI+hDZlf+UUlfHfwaulvReSXtJmi3pJEmfBx4gu4r+ZEo/DlgKfH0aL38/sB04N+V9H5N3RU3WxiL12IWkA8i6hc4B/gxYrNd/SzAf+Je0/nOy7yp2EREnpe8lmi0nNSlvVvryeRAYlDR34ictSb9L9inGd/lUhIO/dctngU+lbp2/aNwQERuAy8mC7bNkXRn3TveFI+Jy4HzgU2RvIE+TfYpYnb4AXQqcBIwB/x34QEQ8MY3X3Qa8j+yL158BpwGrcrYxdz0aSfoN4HbgCxHxjYj4JfBXwGfSLi8Cv5HW35DqXdSngF+R3U77J2n9UxP2+SCwKiJe6kB51gXyl/Jm/UPSYcD5EfEBSX8J/HNE3FRytawH+crfrI9ExPeBZyV9FzgE+Ltya2S9ylf+ZmY15Ct/M7MacvA3M6shB38zsxqqxKieo6OjsWbNmrKrYWZWNS1/KNmVK39Je0q6R9Jgi+1zJK1t9RP9sbGxma2gmfWvzZvhkENg1qzscfPmsmvUE7rV7XMW2Q9AdjTbmH4AcyfZj2jMzDpn6VJ44gnYsSN7XLq07Br1hNzBX9JNaczxByVtlXTyJLsvA25pyHuXpBPT+mXKZv1ZnfYzM+ucjRthZxpPb+fO7LkVuvJfAmyOiKPJgvbyZjulMdGHImJLQ/Jy4GJJy4DDgfPIBvI6qkB9zMx2d/DBMJBC3cBA9tzyBf80yNP+ZJNuAGwA5kuaJ+l6SdekwA7ZsLwvNuaPiLVkX0ScD5weETtSl9A2SfvkqZOZWVO33gqLFrGdQVi0KHtuue/2ORTYFBGvpOdHAI+SDYJ1c0Tcmqa8u5FsEKhdpqNLw9EuAJ6fMBDUHsArmJl1ytAQrF/PbEGsL7syvSNvt88S4MA0tOs8sk8AVwALyUZUBNgBEBEvkIaBBUhT891INrn2y5JGU/r4jEyv5m2MmZlNT5Hgv4psjPKHgBURcS/wDNkbwMTXvgM4VtJeKd8FEfE4cCmvf1dwPHBbzvqYmVkb8nb7LAHOjohzJ6SvAr6Y7vxp7Fi7GvhERHyHbMJt4LW+//HnZ5CNF25mZjMsb/A/iGyC7V1ExC+ADzVJfyTd3jnY7F7/dEfQ6oh4Mmd9zMysDbmCf0QsnHqv3fJcN8m2bcANeepiZmbt88BuZmY15OBvZlZDDv5mZjXk4G9mVkMO/mZmNeTgb2ZWQw7+ZmY15OBvZlZDDv5mZjXUv8G/6LydRfLXLW9V6+021+N4FdHPbY6Inl9GRkaibcPDEQMDEZA9Dg93L3/d8la13m5zPY5XAm1nqXybY5K4qojo7LvJDJD0HLC1nTwjMDIxbR2s60b+uuUts2y3uRp5yy47rz5o81hEjDbbUIngb2ZmndW/ff5mZtaSg7+ZWQ05+JuZ1VDemby6anR0NNasWVN2NczMqkatNnTlyl/SnpLukTTYYvscSWslNX0zGhsbm9kKmllvq+p9/j2sW90+ZwGrosn8vfDaNI53Aqd1qT5mViVLl8ITT8COHdnj0qX9XW4X5A7+km6StFLSg5K2Sjp5kt2XAbc05L1L0olp/TJJVwGr035mZrvauBF27szWd+7MnvdzuV1Q5Mp/CbA5Io4mC9rLm+0kaQ4wFBFbGpKXAxdLWgYcDpwHPAYcVaA+ZtavDj4YBlK4GhjInvdzuV2QK/hLmgvsD1ySkjYA8yUNSbpW0s0Nu+8HvNiYPyLWkn0RcT5wekTsSF1C2yTtk6dOZtbHbr0VFi1iO4OwaFH2vJ/L7YK8d/scCmyKiFfS8yOARyNiM/DhCcH/V8DcxsySFgMLgOcj4qWGTXsAr2Bm1mhoCNavZ7Yg1teg3C7I2+2zBDhQ0lxJ88g+AVzRbMeIeAEYTJ8WkLQAuBE4FXhZ0mhK35dsHIpXc9bJzMymqUjwXwU8ADwErIiIeyfZ/w7gWEl7pXwXRMTjwKW8/l3B8cBtOetjZmZtKBL8r46IJRExHBFfhezqXdKXgMMlXdSw/9XAByPilxFxTER8G7K+/4g4Ju1zBvA3OetjZmZtyNvnfxCwaWJiRDwPnNMk/ZF0e+dgs3v90x1BqyPiyZz1MTOzNuQK/hGxMEee6ybZtg24IU9dzMysfR7Yzcyshhz8zcxqyMHfzKyGHPzNzGrIwd/MrIYc/M3MasjB38yshhz8zcxqyMHfzKyG+jf4F517s0j+uuWtar3d5uocryLKnIe3zOM9lYjo+WVkZCTaNjwcMTAQAdnj8HD38tctb1Xr7TZX53gl0HaW8sotWnYH6h2TxFVFRGffTWaApOeAre3kGYGRiWnrYF038tctb5llu83VyNuJ/HmVVW7RsjtU77GIGG22oRLB38zMOqt/+/zNzKwlB38zsxpy8Dczq6GuBH9Je0q6R9Jgi+1zJK2V1HRymdHR0QC8ePHixUt7S0vduvI/C1gVTaZwBMZn8roTOK3Z9rGxsRmsmpn1tTLv8+9huYO/pJskrZT0oKStkk6eZPdlwC0Nee+SdGJav0zSVcDqtJ+ZWecsXQpPPAE7dmSPS5eWXaOekHcCd4AlwC0RcZqkY4EvALdN3ClNzj4UEVsakpcDn5b0RuBw4JSUflSB+piZ7W7jRti5M1vfuTN7bvmu/CXNBfYHLklJG4D5kt4r6Zr0ieBdadt+wIuN+SNiLSDgfOD0iNiRuoS2SdonT53MzJo6+GAYSKFuYCB7brm7fQ4FNkXEK+n5EcCjEbE6Ij4CnMPr/fe/AuY2Zpa0GFgAbIuIlxo27QG8gplZp9x6KyxaxHYGYdGi7LnlDv5LgAMlzZU0j+wTwBUN2z8FXA0QES8Ag+nTApIWADcCpwIvSxpN6fuS/RT51Zx1MjPb3dAQrF/PbLbD+vXZcysU/FcBDwAPASsi4l5lPgd8KyIeadj/DuBYSXulfBdExOPApWT9/wDH0+Q7AzMz67y8X/guAc6OiHMnpH8cOAF4g6TfiYgvpfSrgU9ExHeAY8Z3Tn3/48/PAC7MWR8zM2tD3uB/ELBpYmJEXAlc2ST9kXR752Cze/3THUGrI+LJnPUxM7M25Ar+EbEwR57rJtm2DbghT13MzKx9HtvHzLqjjr+07eE2O/ibWXfU8Ze2PdxmB38z6446/tK2h9vs4G9m3VHHX9r2cJsd/M2sO+r4S9sebnMl5vA98sgj4+GHHy67GmbWARKUEXbKKrfkstVqg6/8zcxqyMHfzKyGHPzNzGrIwd/MrIb6N/gX/WVdkfx1y1vVervN3T1eRZRVdj+3OSJ6fhkZGYm2DQ9HDAxEQPY4PNy9/HXLW9V6u83dPV4J5MjUgbLLKrfkslvG1Urc6inpOWBrO3lGYGRi2jpY1438dctbZtluczXyFlVW2X3Q5rGIGG22oRLB38zMOqt/+/zNzKwlB38zsxpy8Dczq6GuBH9Je0q6R9Jgi+1zJK2V1HRmsdHR0QC8ePFS12Xz5uCQQ4JZs7LHzZv7u9zOLS1168r/LGBVNJm/FxifxvFO4LRm28fGxmawambW88qaFKWHJ2MpKnfwl3STpJWSHpS0VdLJk+y+DLilIe9dkk5M65dJugpYnfYzM9tVWZOi9PBkLEUVufJfAmyOiKPJgvbyZjtJmgMMRcSWhuTlwMWSlgGHA+cBjwFHFaiPmfWrsiZF6eHJWIrKFfwlzQX2By5JSRuA+ZLeLulLkm6W9B/Ttv2AFxvzR8RasnGmzwdOj4gdqUtom6R98tTJzPpYWZOi9PBkLEU1/YJ1Gg4FNkXEK+n5EcCjEfE4cI6kAeAGYAXwK2BuY2ZJi4EFwPMR8VLDpj2AVzAzazQ0BOvXM1sQ62tQbhfk7fZZAhwoaa6keWSfAK4AkHQKcBtwO0BEvAAMpk8LSFoA3AicCrwsaTSl70v2U+RXC7THzMymoUjwXwU8ADwErIiIewEi4hsRcRK7fnl7B3CspL1SvgvSp4RLef27guPJ3jTMzGyG5e32WQKcHRHnNiZKOg54H1n3ze0Nm64GPhER3wGOGU9Mff/jz88ALsxZHzMza0Pe4H8QsGliYkTcDdzdJP2RdHvnYLN7/dMdQasj4smc9TEzszbkCv4RsTBHnusm2baN7AtiMzPrAo/tY2ZWQw7+ZmY15OBvZlZDDv5mNj1lTmZeRzN8vB38zWx6+niEy540w8fbwd/MpqePR7jsSTN8vB38zWx6+niEy540w8fbwd/MpqePR7jsSTN8vBUx6UxfPeHII4+Mhx9+uOxqmBkgQVlho6yyK9xmtdrgK38zsxrq3+Bf9DapIvnrlreq9Xabu3u7Zlll17HN0xERPb+MjIxE24aHIwYGIiB7HB7uXv665a1qvd3m9tscWdZcyiq7jm1uyN5qqUSfv6TngK3t5BmBkYlp62BdN/LXLW+ZZbvN3ctbVFll17HNDcYiYrTZhkoEfzMz66z+7fM3M7OWHPzNzGoo70xeXTU6Ohpr1qwpuxpmZlVT/n3+kvaUdI+kwUn2mSNpraRd3pTGxsZmvoJm1p96+XbLEnWz2+csYFU0mcN3XGTTOd4JnNa1WplZf/NopE21Ffwl3SRppaQHJW2VdHIb2ZcBtzS81l2STkzrl0m6Km1anfY1MyvOo5E21e6V/xJgc0QcTRagl08nk6Q5wFBEbGlIXg5cLGkZcDhwXkp/DDiqzXqZWa8rq/ulzNFIe7jLadrBX9JcYH/gkpS0AZgvaUjStZJuTvvNk3S9pGtSYAfYD3ix8fUiYi3ZlxHnA6ePdwelx22S9inQLjPrNWV1v5Q5GmkPdzm1c+V/KLApIl5Jz48AHo2IzRHx4Yb93gfcHBEfAU5Jab8C5ja+mKTFwAJgW0S8NKGsPYBXMLP+UVb3y9AQrF/PbLbD+vXZ827p4S6ndoL/EuBASXMlzSP7BHBFk/0WAk+n9fGr+ReAwfTpAUkLgBuBU4GXJb3282NJ+5L9JPnVdhtjZj2sjpPB9HCb2w3+q4AHgIeAFRFxb5P9niF7A5j4+ncAx0raK73OBRHxOHApu353cDxwWxv1MrMqqONkMD3c5mmP7SPpHuDsiNg4IX1f4DPAicD/BK4EvkjWbfO9iLgx7XcE8ImIOHOKclYBF0bEk+NpnszFrH94QpbuFt1qQzu/8D0I2DQxMSKeB86ZkPyhJvs9km7vHGx1r3+6K2h1Y+A3M7POm3bwj4iFU+815WtcN8X2bcANRcsxM7PJeWA3szrp4fvObYIZPlcO/mZ10sP3ndsEM3yuHPzN6qSH7zu3CWb4XDn4m9VJD993bhPM8Lly8Derkx6+79wmmOFzVYk5fH2fv1ln1fGe96q2uWC9y5/Mxcwa+K6baujj8+Tgb1YG33VTDX18nhz8zcrgu26qoY/Pk4O/la+sj9ZlfqSv4l03ZR6vOk4EM9MioueXkZGRaNtTT0UMD0cMDmaPTz3Vvfx1y1s0//BwxMBABGSPw8O9X27RslPeVynhXEXW5LYVPV4ll52r3CLnqWjZHcgbk8TV0gP7dJZcwb/oH0uR/HXLWzT/4GCWb3wZHOz9couWnVQqEBY9XiWXXWIA7sngX4lbPSU9B2xtJ88IjExMWwfrupG/bnmL5j8UDtmjYaa3X8Mrj8H6Xi63aNlFlFVu0eNV1bIrbiwiRpttqETwNzOzzvIXvmZmNeTgb2ZWQ+3M5FWa0dHRWLNmTdnVMDOrmt4a3kHSnpLukTTYYvscSWslzQIYGxvrbgXNrPPq+HuOHh4eoqxun7OAVdFiLt/IpnO8Ezitq7Uy63dlBqOyhkooc4iGImXP8Lma1t0+km4ie6N4K/Am4M8i4rbchUr3AWdExBZJdwH/NSK+Leky4A0R8XFJS4DPRsS7PaqnWYccckgWhHbuzH6xumgRrO/SHZOzZmVBcNzgIGzf3r/lFi27M+eqcLfPEmBzRBwNLAOWt1uD12oizQGGImJLSloOXCxpGXA4cF5Kfww4Km85Zj2trCvwMseqKWuohCLlFj1PRcoueyYvSXOB/YFLUtIGYL6kIUnXSro57TdP0vWSrkmBvGkasB/w4vjrR8Rasnen84HTx7uC0uM2Sft0pKVmvaSsrogyx6opayKZIuUWPU9Fyp7pczXZz39Tl9CRwH0Nz08Abm54fnN6PBNYmtZXTpI2H9jSkH8x8CRwf5Oyx4DZuYZ3MJtKB8bIya1DQyW0rUNj1RRRZIydro9nVOaQFp05Vy1j+3S6fZYAB0qaK2ke2SeAK5rstxB4Oq3vaJUWES8Ag+n1FgA3AqcCL0t67WfIkvYl+2nyq9Ooo1n7yvwisKwr8KEhWL+e2WzP+o+HhrpTblF1/KQ0w+dqusF/FfAA8BCwIiLubbLfM2TBvvF1m6UB3AG8K73uBRHxOHApu36XcDyQ+0tlsymV2f/tuXTbU9a56uPzNN3gf3VELImI4Yj4KmRX5pK+BBwu6SKyQP5+SSuA8SPULA3gauDfR8QxEfFtyPr+I+KYhn3OAP6mUOusGuo4VnsVr8DrOP9BFc/TNE15q6ekZ4ADI2JnRwuWzgKujyb3+qc7gk6PiBvAE7j3vbJuP9y8GZYuZfuGjcwaPji7quvyP3elJjPv0HnKVXYHzlWJk6j35ATulRjV08G/C9I/Fxs3ZldV3QyEZd6HTXkBOHfZHThXucrt0HkqKxA6+O/KA7tZpo5fflZVHb/8tI5z8O81dfzxT1W/VKvbuarqebKm3O3Ta8rq/y7zZ/9Jpfq/oSPHrKy+9zp2gdSxzbjbJwdf1XWn3KLKvAPF58oqrH+Df9GgULd+1aK3tJUVhOv4XUUf335o3dO/wb9oUChyVVckEBa5qqvjcL3+rsIsl/7t8y96W1qRftU69gV34DbAyrW5A/nrlrfMsuvYZmrZ51/0I3mRq7qyrkbrOFyvr77Ncunf4F80KBTpV63iuOVFlRWE3f9tlkv/Bv8yg0IVxy0vykHYrFL6t88/qWpfXRXzllm221yNvGWWXcc2U/WxfSQ9B2wtux5mZhUzFhGjzTZUIvibmVln9W+fv5mZteTgb2ZWQw7+ZmY1VErwl7SnpHskDbbYPkfSWkmzAEZHRwPw4sWLFy/tLS2VdeV/FrCq2RSOABGxDbgTOA1gbGysi1Uz62Nljv9UVtlVbfNM1zsiplyAm4CVwINkt1yePJ18k7zefcBb0vpdwIlp/TLgqrS+BLg9IhgZGQmzvvLUUxHDwxGDg9njU091p9zh4YiBgQjIHoeHu1NumWUXKbfoeSpSdmeOV+s4PNnGeD1YbwA+m9aPBR6cTr4WrzUH+EnD83cCdwPLgNuAwZQ+CDwXDv42U8oKwBHlBcLBwazM8WVwsDvlFi27yLkqUm7R81Sk7M6cq5axeMpuH0lzgf2BS1LSBmC+pPdKukbSSknvkjRP0vUpbVnKu1sasB/wYsMnj7Vkv0I7Hzg9UldQetwmaZ+pP79YpdVxLoCyBuErMv5T0fNUpOwi56pIuUXPU5GyZ3qsrsneGSK7Aj8SuK/h+QnAzQ3P5wPXAmcCS1PayvTYLG0+sKUh/2LgSeD+JmWPAbN95d8Fvgru7lVwWV0RKe+r5Mhb9DwVKbsDnxoq1+YieV+Xv9sH+DDwDDAXmAfcC/xew/bLgSOAi4DDUtrX0uNuaWn96fR6C4AfAG8Hvg2MNuyzL/BE1K3bp459wWUF4Q79Y1cuCEf6z29Xh85TrrLLanNnAnC+sjuQNwoG/yvT8ihZl8+fpHQBnwNOiNev8t+T1r/eKi2tXwucAtzP61/2vrPx6h/4Q+DyyBv8iwbRugXhov/YRY5XFa+Ci9Y7KSsIlxWAc5fdgSBcYgCubPC/Bzi4Sfq5wDrgS8A56VPB/wJWAMvSPrulpfQjgK9MUe4q4F9H3uBf9A+1rIBU1pdiPXC8KnUVHFHpIFzHq2AH//aD/zPAwFT7tbuQ3es/2GLbHOAD489zBf+i/5hl3SFQVt6i/9hlBcKyyo2odBCuYyCsY5tjkhg85d0+EbEwInZOtV+7IuK6mORHXhFxQ6ECin5TXtYdAmVNH1l0MhbPXubZy6xS+ndsn6L/mEXyFwlIVZw+Ejx7mYOwVUwlxvOv3ExemzfD0qVs37CRWcMHZwGpzcBQqXI7lL+Kecssu4p5yyy7jm2m6jN5VS741zRvmWW7zdXIW2bZdWwzkwT//u32MTOzlhz8zcxqyMHfzKyGHPzNzGrIwd/MrIYc/M3MasjB38yshhz8zcxqyMHfzKyGKvELX0nPkU0cb2Zm0zcWEaPNNlQi+JuZWWe528fMrIYc/M3MasjB38yshkoJ/pL2lHSPpMEW2+dIWitpFsDo6GgAXrz0z7J5c3DIIcGsWdnj5s39XW5V21y0zmWWnS2tTTbH40wtwEeBP59in+WkSd9zzeFbZDLzovnrlreq9S6zzR2Yb7lS5ZZZdlnzYpdddqZ1jJ1sY7weiG8CVgIPkt1yefJ08k3yevcBb0nrdwEnpvXLgKvS+hLg9sgb/Kt60qqYt6r1LrPNRSaer2K5ZZZdpNwiecsuO1M4+G8APpvWjwUenE6+Fq81B/hJw/N3AncDy4DbgMGUPgg8F3mDf1VPWhXzVrXeZba5igG86PGqYpurelH0upaxeMo+f0lzgf2BS1LSBmB+2jZP0sOS3pPWr5d0jaRlDdt3SQP2A15s6HZaSzbV2PnA6RGxI6XvALZJ2meqOjZVdDLzIvnrlrfMsqva5iITz1ex3DLLLlJukbxllz2Vyd4ZIrsCPxK4r+H5CcDNaf3TwCeB9wBnAktT+sr02CxtPrCl4fUWA08C9zcpewyYXaTP/1WK9efmyl+3vFWtd5ltTqDtLJUut8yyc5Xbgbwll90ytk/5C19JHya76v8dsq6YO1LA3wvYF5ibgvRi4FsR8X1JX4uIMyRdNDEtvebTwNvSG8E/AKcBVwKXR8SatM++wL0RscgTuFcjb5llu83VyFtm2XVsMwUncF8CrAIeAB4CVkTEvcBxwDuAM4CPAD8CFk543WeapEH2BvKu9LoXRMTjwKVkd/iMO57sOwAzM+uwWdPYZwlwdkSc25gYERcDSPpTsiv/u4AvSjoZGO+cWtUkDeBq4BMRcUzD660FjmnY5wzgwrZaY2Zm0zKd4H8QsKnVxoj4csPTD03Y9ouJaSn9EUl3SRqM9AVvI0lzgNUR8eQ06mdmZm2aMvhHxMKp9skjIq6bZNs24IaZKNfMzDy2j5lZLTn4m5nVkIO/mVkNOfibmdWQg7+ZWQ05+JuZ1ZCDv5lZDTn4m5nVkIO/mVkNOfibmdWQg7+ZWQ05+JuZ1dCUk7n0AknPkU0cn8d+ZENO2/T4eLXPx6w9Pl7tKXK8xiJitNmGSgT/IiQ9HBFHll2PqvDxap+PWXt8vNozU8fL3T5mZjXk4G9mVkN1CP7/o+wKVIyPV/t8zNrj49WeGTlefd/nb2Zmu6vDlb+ZmU3g4G9mVkN9HfwljUraKOmfJF1Ydn16naQtkv5R0vclPVx2fXqNpOsk/VTSYw1pvyXp25I2pcf5Zdaxl7Q4Xv9F0g/T39j3Jb27zDr2EklvlnSXpA2S1kv685Q+I39jfRv8JQ0CVwMnAcPAH0saLrdWlXB8RBzm+7Cb+jIw8QczFwJ3RsTbgDvTc8t8md2PF8AV6W/ssIi4vct16mXbgQsiYhh4B/DRFLNm5G+sb4M/cDTwTxGxOSK2AV8HTi25TlZhEbEW+NmE5FOB69P69cB7u1mnXtbieFkLEfHjiHgkrb8EPA4cwAz9jfVz8D8AeLrh+TMpzVoL4A5J6ySdXXZlKuJNEfHjtP4T4E1lVqYiPibpB6lbyN1kTUh6C3A48AAz9DfWz8Hf2ndsRBxB1lX2UUnvLLtCVRLZfdO+d3pyK4CDgMOAHwOXl1qbHiRpb+DvgPMi4l8at3Xyb6yfg/8PgTc3PF+Y0qyFiPhhevwp8PdkXWc2uWclLQBIjz8tuT49LSKejYgdEbETuAb/je1C0myywH9jRKxKyTPyN9bPwf8h4G2S3ippDnA68I2S69SzJM2TtM/4OvAu4LHJcxnZ39QH0/oHgVtKrEvPGw9iyR/gv7HXSBJwLfB4RHyhYdOM/I319S98021kfw0MAtdFxGfKrVHvkjREdrUPMAv4mo/XriTdBBxHNsTus8ByYDXwv4EDyYYd/6OI8JectDxex5F1+QSwBfgPDf3ZtSbpWOC7wD8CO1PyX5L1+3f8b6yvg7+ZmTXXz90+ZmbWgoO/mVkNOfibmdWQg7+ZWQ05+JuZ1ZCDv5lZDTn4m5nVkIO/WQ5p3PUT0/plkq4qu05m7ZhVdgXMKmo58GlJbyQbffGUkutj1hb/wtcsJ0n3AHsDx6Xx180qw90+ZjlIWgwsALY58FsVOfibtSmNTHkj2QxLL0tqNlWhWU9z8Ddrg6S9gFVkc60+DlxK1v9vVinu8zczqyFf+ZuZ1ZCDv5lZDTn4m5nVkIO/mVkNOfibmdWQg7+ZWQ05+JuZ1dD/B3MX3poNwDFQAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(len(states), 1)\n",
    "axs[0].set_title(f\"Initial Condition $x_0 = {x_0}$\")\n",
    "for ind, state_no in enumerate(states):\n",
    "    ax = axs[ind]\n",
    "    q = calculate_q_at_state_no(state_no, p_initial, transition_mat)\n",
    "    plot_stem(ax)\n",
    "    if state_no != states[-1]:\n",
    "        axs[ind].set_xticks([])\n",
    "    else:\n",
    "        axs[ind].set_xticks(range(0, 21, 5))\n",
    "        axs[ind].set_xlabel(\"$x$\")\n",
    "\n",
    "pml.savefig(f\"randomWalk{x_0}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d5597ff1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.99999994 0.047619047619047616 [0.04665894 0.04683788 0.04656532 0.04708623 0.04656531 0.04738191\n",
      " 0.04665894 0.04769868 0.04683788 0.04800836 0.04708624 0.04828346\n",
      " 0.04738192 0.04849952 0.04769868 0.04863734 0.04800837 0.04868469\n",
      " 0.04828347 0.04863734 0.04849952]\n"
     ]
    }
   ],
   "source": [
    "print(sum(q), 1 / 21, q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47368c19",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:py3713_2]",
   "language": "python",
   "name": "conda-env-py3713_2-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.13"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}