parte_9.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7660b2c7",
   "metadata": {},
   "source": [
    "# 9 - Markov Chain Monte Carlo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2cea1405",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from scipy import stats\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.gridspec import GridSpec\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "import networkx as nx\n",
    "# from causalgraphicalmodels import CausalGraphicalModel\n",
    "\n",
    "import arviz as az\n",
    "# ArviZ ships with style sheets!\n",
    "# https://python.arviz.org/en/stable/examples/styles.html#example-styles\n",
    "az.style.use(\"arviz-darkgrid\")\n",
    "\n",
    "import xarray as xr\n",
    "\n",
    "import stan\n",
    "import nest_asyncio\n",
    "\n",
    "plt.style.use('default')\n",
    "plt.rcParams['axes.facecolor'] = 'lightgray'\n",
    "\n",
    "# To DAG's\n",
    "import daft\n",
    "from causalgraphicalmodels import CausalGraphicalModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "449b2877",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add fonts to matplotlib to run xkcd\n",
    "\n",
    "from matplotlib import font_manager\n",
    "\n",
    "font_dirs = [\"fonts/\"]  # The path to the custom font file.\n",
    "font_files = font_manager.findSystemFonts(fontpaths=font_dirs)\n",
    "\n",
    "for font_file in font_files:\n",
    "    font_manager.fontManager.addfont(font_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2774f695",
   "metadata": {},
   "outputs": [],
   "source": [
    "# To make plots like drawing \n",
    "# plt.xkcd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6f4405d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# To running the stan in jupyter notebook\n",
    "nest_asyncio.apply()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c807993a",
   "metadata": {},
   "source": [
    "### R Code 9.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0a0ff53a",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_weeks = 10000\n",
    "\n",
    "positions = np.zeros(num_weeks, dtype=int)\n",
    "\n",
    "current = 10\n",
    "\n",
    "for i in range(num_weeks):\n",
    "    # Record current position\n",
    "    positions[i] = current\n",
    "    \n",
    "    # F/lip coin to generate proposal\n",
    "    proposal = current + np.random.choice([-1, 1], size=1)\n",
    "    \n",
    "    # Now make sure he loops around the archipelago\n",
    "    if proposal < 1: \n",
    "        proposal = 10\n",
    "    \n",
    "    if proposal > 10:\n",
    "        proposal = 1\n",
    "    \n",
    "    # Move?\n",
    "    prob_move = proposal / current\n",
    "    \n",
    "    if np.random.uniform(0, 1) < prob_move:\n",
    "        current = proposal\n",
    "    else:\n",
    "        current = current"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e922932",
   "metadata": {},
   "source": [
    "### R Code 8.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "793936d3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+QAAAGDCAYAAACiM9raAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAAAikUlEQVR4nO3dfbRleVkf+O9TdaVLIEApPXdsQEFDpCnuqKSiqJGwoDODAYdZWUZBCYbBMMkYAV9icHqMZjI9mhWDwhiZ1SKio4PVgywlzCzRRhl8g9gtajVc30Z5by4XKVrA1IWqeuaPe3osLreqblWfe377nvv5rNWrztl7n72fffZv//b59n651d0BAAAAFuvI6AIAAADgMBLIAQAAYACBHAAAAAYQyAEAAGAAgRwAAAAGEMgBAABgAIEcAOakqj63qj5WVUdH13I1quqdVXXTwOW/vaqedIVpDuR3CwCXI5ADwFXYGV6r6plVdaaq/k53v7u7H9jd5/dhua+qqq6qZ+wY/kOz4f9o3stclO4+0d1vusI0n/LdVtWbquqbF1IgAOwTgRwArlFVfVOSf5/kad39/yxgkX+U5DkXLX8lydcl+X+vZWazzwMAgwjkAHANquq/S/LvkvxX3f2bs2GPnJ2tXpm9f1NV/euq+o2q+mhV/VJVPfSieTynqt5VVX9eVd+zh0vH/0OSv11Vx2fvn5rk95N84KJ5fkFV/cpsnh+qqp+pqodcNP6dVfUvqur3k3x8Zyivqhur6s+q6lmz9/+4qv6kqj5cVa+rqhtmw19eVT+447O/UFXfvst3ddlpL17vqvrSqrqjqv6iqjaq6iU7v9uquiXJVyX5kdll7D9S236oqj44++zpqnrcZb5LABhOIAeAq/dPk/xPSZ7S3XdcYdpvSPLcJP9Zkvsl+c4kqarHJvnRJN+Y5HOSPDjJw64wr7NJfiHJM2fvn5Pkp3ZMU0m+P8kNSW5M8ogk37djmmcleVqSh3T3uf//g1WPT/KGJN/a3a+uqifP5vV1sxrfleRnZ5O/OsnXV1XNPns8yX950fiLXc20L03y0u5+UJIvSHLbzgm6++Ykv5bkn80uY/9ns/k9McnfyPZ3+XVJ/nyX+QPAZAjkAHD1/m6StyQ5vYdpf6K7/6i7/1O2w+UXz4Z/bZL/0N2/3t2fSPIvk/Qe5vdTSZ4zO+v9d5L8/MUju/tPuvuXu3uruzeTvGQ23cVe1t3vmdV0r69K8rokz+nu18+GfWOSV3b373T3VpLvTvLlVfXIbAfinn3u3vX5re5+/y41X820n0zy16vqod39se5+y+W+jB2f+2tJHpOkunu9u+/e42cBYAiBHACu3j/N9pnYV9x71vcyPnDR679M8sDZ6xuSvOfeEd39l9nDGd3u/vUk1ye5Ocnrd4TqVNVqVf1sVb2vqv4iyU8neeiO2bwnn+6fJPnNHQ9XuyHbZ8XvXfbHZjU+rLs722e4nzUb/Q1JfuYSNe952iTPy/Z3+wdV9dtV9fRLTLdzGb+S5EeyfU//B6vq1qp60F4+CwCjCOQAcPU2kjwl22d8f/Qa53F3koff+6aqPjPJZ+/xsz+d5Dvy6ZerJ8n/ku2z0Wuzy76fne3L2C+225n4f5Lkc6vqhy4a9v4kn3dRjQ+Y1fi+2aBXJ/naqvq8JF+W5OcuU/Oepu3uP+7uZ2X7Ev9/k+Q1s+V+2qS7fPZl3f03kzw226H+n1+mHgAYTiAHgGswu9z6KUmeuiPE7tVrknxNVX1FVd0v2/d5X+ls+71elu3L5t+8y7i/luRjSe6pqodl76H0o9l+SNwTq+oHZsNeneS5VfXFVXVdtsP+W7v7nUnS3W9L8qEkr0jyhu7+yKVmvtdpq+rZVXV9d19Icu80F3aZdCPJ51/0ub9VVV9WVZ+R5OPZvt9+t88BwGQI5ABwjbr73UmenO0zv99/lZ99e5Jvzfal3HdnO0R/MMnWHj774e5+4+xS8J3+VZLHJ7knyf+V5LVXUdNHsh30v7qq/nV3357ke7J9NvvubD9k7Zk7PvZ/JLlp9u+V7GXapyZ5e1V9LNsPeHvmzsvyZ16a7e/9TFW9LMmDkvxYkjPZvsz+z5P82z3UBADD1O7HcgBgkarqgdk+I/zo7v6zweUAAAvgDDkADFJVX1NV95/dI/2D2X5q+zvHVgUALIpADgDjPCPbD057f5JHZ/vybJeuAcAh4ZJ1AAAAGMAZcgAAABhAIAcAAIABVkYXsBfHjx/vG264YXQZAAAAcFXe8Y53fKi7r99t3IEI5DfccENOnTo1ugwAAAC4Kmtra++61DiXrAMAAMAAAjkAAAAMIJADAADAAAI5AAAADCCQAwAAwAACOQAAAAwgkAMAAMAAAjkAAAAMIJADAADAAPsWyKvqlVX1waq666Jhn1VVv1xVfzz79/h+LR8AAACmbD/PkL8qyVN3DHtxkjd296OTvHH2fmmsb27ltrvuyfrm1jVPM495LKrWg1bvQVtnAODwmtJvrEXOh+U3pXY7BSv7NePufnNVPXLH4GckedLs9U8meVOSf7FfNSzS+uZWbr59I+fOd1aOVm65aTU3Xn/dVU0zj3ksqtaDVu9BW2cA4PCa0m+sRc6H5TeldjsVi76HfLW77569/kCS1UtNWFXPr6o7quqOM2fOLKa6++D0xtmcO9+5kOTchc7pjbNXPc085rGoWg9avQdtnQGAw2tKv7EWOR+W35Ta7VQMe6hbd3eSvsz4W7v7ZHefPH58+rear60ey8rRypFKVo5U1laPXfU085jHomo9aPUetHUGAA6vKf3GWuR8WH5TardTUdu5eJ9mvn3J+uu7+3Gz93+Y5EndfXdVfU6SN3X3F15pPidOnOhTp07tW53zsr65ldMbZ7O2euySl01caZp5zGNRtR60eg/aOgMAh9eUfmMtcj4svym120VZW1u7s7tP7jZu0YH83yb58+7+gap6cZLP6u7vutJ8DkogBwAAgItdLpDv5589e3WS30ryhVX13qp6XpIfSPJ3q+qPk9w0ew8AAACHzn4+Zf1Zlxj1lP1aJgAAABwUwx7qBgAAAIeZQA4AAAADCOQAAAAwgEAOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQA4AAAADCOQAAAAwgEAOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQA4AAAADCOQAAAAwgEAOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQA4AAAADCOQAAAAwgEAOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQM5w65tbue2ue7K+uTW6lIWZ1zofxu8OgOnYy3FoHseqZTveLWp9pvR7Y1FtheVwmNrCyugCONzWN7dy8+0bOXe+s3K0cstNq7nx+utGl7Wv5rXOh/G7A2A69nIcmsexatmOd4tanyn93lhUW2E5HLa24Aw5Q53eOJtz5zsXkpy70Dm9cXZ0SftuXut8GL87AKZjL8eheRyrlu14t6j1mdLvjUW1FZbDYWsLAjlDra0ey8rRypFKVo5U1laPjS5p381rnQ/jdwfAdOzlODSPY9WyHe8WtT5T+r2xqLbCcjhsbaG6e3QNV3TixIk+derU6DLYJ+ubWzm9cTZrq8eW+nKUi81rnQ/jdwfAdOzlODSPY9WyHe8WtT5T+r2xqLbCcli2trC2tnZnd5/cbZxADgAAAPvkcoHcJesAAAAwgEAOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQA4AAAADCOQAAAAwgEAOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQA4AAAADCOQAAAAwgEAOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQA4AAAADCOQAAAAwgEAOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQA4AAAADCOQAAAAwgEAOAAAAAwwJ5FX1bVX19qq6q6peXVXHRtQBAAAAoyw8kFfVw5K8IMnJ7n5ckqNJnrnoOmDq1je3cttd92R9c2voPADYm2Xrcxe1PvNazpSOm8vWFqZkStsZ5mFl4HI/s6o+meT+Sd4/qA6YpPXNrdx8+0bOne+sHK3cctNqbrz+uoXPA4C9WbY+d1HrM6/lTOm4uWxtYUqmtJ1hXhZ+hry735fkB5O8O8ndSe7p7l/aOV1VPb+q7qiqO86cObPoMmGo0xtnc+5850KScxc6pzfODpkHAHuzbH3uotZnXsuZ0nFz2drClExpO8O8jLhk/XiSZyR5VJIbkjygqp69c7ruvrW7T3b3yePHjy+6TBhqbfVYVo5WjlSycqSytnr1j1mYxzwA2Jtl63MXtT7zWs6UjpvL1hamZErbGealunuxC6z6B0me2t3Pm71/TpIndPd/f6nPnDhxok+dOrWoEmES1je3cnrjbNZWj13zpVTzmAcAe7Nsfe6i1mdey5nScXPZ2sKUTGk7w16tra3d2d0ndxs3IpB/WZJXJvlbSf5TklcluaO7/9dLfUYgBwAA4CC6XCAfcQ/5W5O8JsnvJDk9q+HWRdcBAAAAIw15ynp3f2+S7x2xbAAAAJiChZ8hBwAAAARyAAAAGEIgBwAAgAEEcgAAABhAIAcAAIABBHIAAAAYQCAHAACAAQRyAAAAGEAgBwAAgAEEcgAAABhAIAcAAIABBHIAAAAYQCAHAACAAQRyAAAAGEAgBwAAgAEEcgAAABhAIAcAAIABBHIAAAAYQCAHAACAAQRyAAAAGEAgBwAAgAEEcgAAABhAIAcAAIABBHIAAAAYQCAHAACAAQRyOOTWN7dy2133ZH1za3QpAJ9iSv3TvGqZyjpNpY5FWrZtyO4WtZ21A+ZlZXQBwDjrm1u5+faNnDvfWTlaueWm1dx4/XWjywKYVP80r1qmsk5TqWORlm0bsrtFbWftgHlyhhwOsdMbZ3PufOdCknMXOqc3zo4uCSDJtPqnedUylXWaSh2LtGzbkN0tajtrB8yTQA6H2NrqsawcrRypZOVIZW312OiSAJJMq3+aVy1TWaep1LFIy7YN2d2itrN2wDxVd4+u4YpOnDjRp06dGl0GLKX1za2c3jibtdVjLrcCJmVK/dO8apnKOk2ljkVatm3I7ha1nbUDrsba2tqd3X1yt3ECOQAAAOyTywVyl6wDAADAAAI5AAAADCCQAwAAwAACOQAAAAwgkAMAAMAAAjkAAAAMIJADAADAAAI5AAAADCCQAwAAwAACOQAAAAwgkAMAAMAAAjkAAAAMIJADAADAAAI5AAAADCCQAwAAwAACOQAAAAwgkAMAAMAAAjkAAAAMIJADAADAAAI5AAAADCCQAwAAwAACOQAAAAwgkAMAAMAAAjkAAAAMIJADAADAAAI5AAAADDAkkFfVQ6rqNVX1B1W1XlVfPqIOAAAAGGXUGfKXJvnF7n5Mki9Ksj6oDmAP1je3cttd92R9c2voPIBpmNf+PJW+ZUr905S+W+C+sR+yFyuLXmBVPTjJE5P8oyTp7k8k+cSi6wD2Zn1zKzffvpFz5zsrRyu33LSaG6+/buHzAKZhXvvzVPqWKfVPU/pugfvGfshejThD/qgkm0l+oqreVlWvqKoH7Jyoqp5fVXdU1R1nzpxZfJVAkuT0xtmcO9+5kOTchc7pjbND5gFMw7z256n0LVPqn6b03QL3jf2QvRoRyFeSPD7Jy7v7S5J8PMmLd07U3bd298nuPnn8+PFF1wjMrK0ey8rRypFKVo5U1laPDZkHMA3z2p+n0rdMqX+a0ncL3Df2Q/aqunuxC6z6z5O8pbsfOXv/VUle3N1Pu9RnTpw40adOnVpQhcBO65tbOb1xNmurx675cqt5zAOYhnntz1PpW6bUP03puwXuG/sh91pbW7uzu0/uNm7hgTxJqurXknxzd/9hVX1fkgd09z+/1PQCOQAAAAfR5QL5wh/qNvOtSX6mqu6X5E+TPHdQHQAAADDEngJ5VVXvOJVeVdd19zU9w7+7fzfJrv+HAAAAAA6DvT7U7ccvflNVD0zyf8+/HAAAADgc9hrI31tVP5okVXU8yS8l+el9qwoAAACW3J4CeXf/yyQfq6r/Ldth/N9190/sa2UAAACwxC57D3lV/f2L3r41yfck+Y9Juqr+fne/dj+LAwAAgGV1pYe6fc2O929L8hmz4Z1EIAcAAIBrcNlA3t3+HBkAAADsgz3dQ15Vf6Oq3lhVd83e/xdV9T/ub2kAAACwvPb6lPUfS/LdST6ZJN39+0meuV9FAQAAwLLbayC/f3f/xx3Dzs27GAAAADgs9hrIP1RVX5DtB7mlqr42yd37VhUAAAAsuSs9Zf1e35Lk1iSPqar3JfmzJM/et6oAAABgye0pkHf3nya5qaoekORId390f8sCAACA5bbXp6yvVtWPJ3lNd3+0qh5bVc/b59oAAABgae31HvJXJXlDkhtm7/8oyYv2oR4AAAA4FPYayB/a3bcluZAk3X0uyfl9qwoAAACW3F4D+cer6rPzV09Zf0KSe/atKgAAAFhyl32oW1W9KMlvJvmuJL+Q5POr6jeSXJ/kH+x7dQAAALCkrvSU9Ycn+eEkj0nyB0l+Ocmbk7y6uz+0v6UBAADA8rpsIO/u70ySqrpfkpNJviLJk5J8d1V9pLsfu+8VAgAAwBLa098hT/KZSR6U5MGz/96f5PR+FQUAAADL7kr3kN+a5ESSjyZ5a7bvJ39Jd59ZQG0AAACwtK70lPXPTXJdkg8keV+S9yb5yD7XBAAAAEvvSveQP7WqKttnyb8iyXckeVxVfTjJb3X39y6gRgAAAFg6V7yHvLs7yV1V9ZFs/+3xe5I8PcmXJhHIAQAA4Bpc6R7yF2T7zPhXJPlktu8h/80kr4yHugEAAMA1u9IZ8kcm+T+TfFt3373/5QAH1frmVk5vnM3a6rHceP11w+YxpeXAok1pP5zSfjalWuZh2dYHDjP7M1e6h/zbF1UIcHCtb27l5ts3cu58Z+Vo5ZabVq/6oDKPeUxpObBoU9oPp7SfTamWeVi29YHDzP5McuWnrANc0emNszl3vnMhybkLndMbZ4fMY0rLgUWb0n44pf1sSrXMw7KtDxxm9mcSgRyYg7XVY1k5WjlSycqRytrqsSHzmNJyYNGmtB9OaT+bUi3zsGzrA4eZ/Zkkqe2HqE/biRMn+tSpU6PLAC5jSveuTmU5sGhT2g+ntJ9NqZZ5WLb1gcPM/nw4rK2t3dndJ3cbJ5ADAADAPrlcIHfJOgAAAAwgkAMAAMAAAjkAAAAMIJADAADAAAI5AAAADCCQAwAAwAACOQAAAAwgkAMAAMAAAjkAAAAMIJADAADAAAI5AAAADCCQAwAAwAACOQAAAAwgkAMAAMAAAjkAAAAMIJADAADAAAI5AAAADCCQAwAAwAACOQAAAAwgkAMAAMAAAjkAAAAMIJADAADAAAI5AAAADCCQAwAAwAACOQAAAAwgkAMAAMAAwwJ5VR2tqrdV1etH1QAAAACjjDxD/sIk6wOXDxxA65tbue2ue7K+uXVN4+e1HMbay/aZ0jacR73aNsDhpN9ebisjFlpVD0/ytCS3JPn2ETUAB8/65lZuvn0j5853Vo5WbrlpNTdef92ex89rOYy1l+0zpW04j3q1bYDDSb+9/EadIf/hJN+V5MKlJqiq51fVHVV1x5kzZxZWGDBdpzfO5tz5zoUk5y50Tm+cvarx81oOY+1l+0xpG86jXm0b4HDSby+/hQfyqnp6kg92952Xm667b+3uk9198vjx4wuqDpiytdVjWTlaOVLJypHK2uqxqxo/r+Uw1l62z5S24Tzq1bYBDif99vKr7l7sAqu+P8k/THIuybEkD0ry2u5+9qU+c+LEiT516tSCKgSmbH1zK6c3zmZt9diul2xdafy8lsNYe9k+U9qG86hX2wY4nPTbB9/a2tqd3X1yt3ELD+SfsvCqJyX5zu5++uWmE8gBAAA4iC4XyP0dcgAAABhgyFPW79Xdb0ryppE1AAAAwAjOkAMAAMAAAjkAAAAMIJADAADAAAI5AAAADCCQAwAAwAACOQAAAAwgkAMAAMAAAjkAAAAMIJADAADAAAI5AAAADCCQAwAAwAACOQAAAAwgkAMAAMAAAjkAAAAMIJADAADAAAI5AAAADCCQAwAAwAACOQAAAAwgkAMAAMAAAjkAAAAMIJADAADAAAI5AAAADCCQAwAAwAACOQAAAAwgkAMAAMAAAjnANVrf3Mptd92T9c2taxq/yOXMq5YrWdQ6z8s8lrOoWgFgNwft2MunWhldAMBBtL65lZtv38i5852Vo5VbblrNjddft+fxi1zOvGq5r7Uuej6LWM6iagWA3Ry0Yy+fzhlygGtweuNszp3vXEhy7kLn9MbZqxq/yOXMq5b7Wuui57OI5SyqVgDYzUE79vLpBHKAa7C2eiwrRytHKlk5UllbPXZV4xe5nHnVcl9rXfR8FrGcRdUKALs5aMdePl119+garujEiRN96tSp0WUAfIr1za2c3jibtdVju17WdaXxi1zOvGq5r7Uuej6LWM6iagWA3Ry0Y+9htLa2dmd3n9xtnEAOAAAA++Rygdwl6wAAADCAQA4AAAADCOQAAAAwgEAOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQA4AAAADCOQAAAAwgEAOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQA4AAAADCOQAAAAwgEAOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQA4AAAADCOQAAAAwgEAOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQA4AAAADLDyQV9UjqupXq+odVfX2qnrhomsAAACA0UacIT+X5Du6+7FJnpDkW6rqsQPqABhufXMrt911T9Y3t0aXsqda5lHvlNZ5Lw5avQCw02E8fh8UK4teYHffneTu2euPVtV6kocleceiawEYaX1zKzffvpFz5zsrRyu33LSaG6+/brK1zKPeKa3zXhy0egFgp8N4/D5Iht5DXlWPTPIlSd66y7jnV9UdVXXHmTNnFl4bwH47vXE25853LiQ5d6FzeuPspGuZR71TWue9OGj1AsBOh/H4fZAMC+RV9cAkP5fkRd39FzvHd/et3X2yu08eP3588QUC7LO11WNZOVo5UsnKkcra6rFJ1zKPeqe0zntx0OoFgJ0O4/H7IKnuXvxCqz4jyeuTvKG7X3Kl6U+cONGnTp3a/8IAFmx9cyunN85mbfXY8Eu/9lLLPOqd0jrvxUGrFwB2OozH7ylZW1u7s7tP7jZu4YG8qirJTyb5cHe/aC+fEcgBAAA4iC4XyEdcsv6VSf5hkidX1e/O/vt7A+oAAACAYUY8Zf3Xk9SilwsAAABTMvQp6wAAAHBYCeQAAAAwgEAOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQA4AAAADCOQAAAAwgEAOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQA4AAAADCOQAAAAwgEAOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQA4AAAADCOQAAAAwgEAOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQA4AAAADCORLan1zK7fddU/WN7dGlwIAABwCMsjVWxldAPO3vrmVm2/fyLnznZWjlVtuWs2N1183uiwAAGBJySDXxhnyJXR642zOne9cSHLuQuf0xtnRJQEAAEtMBrk2AvkSWls9lpWjlSOVrByprK0eG10SAACwxGSQa+OS9SV04/XX5ZabVnN642zWVo+5VAQAANhXMsi1EciX1I3XX2cnAAAAFkYGuXouWQcAAIABBHIAAAAYQCAHAACAAQRyAAAAGEAgBwAAgAEEcgAAABhAIAcAAIABBHIAAAAYQCAHAACAAQRyAAAAGEAgBwAAgAGqu0fXcEVVtZnkXaPr2KOHJvnQ6CLgKmm3HETaLQeVtstBpN1yEE2l3X5ed1+/24gDEcgPkqq6o7tPjq4DroZ2y0Gk3XJQabscRNotB9FBaLcuWQcAAIABBHIAAAAYQCCfv1tHFwDXQLvlINJuOai0XQ4i7ZaDaPLt1j3kAAAAMIAz5AAAADCAQD4nVfXUqvrDqvqTqnrx6HrgUqrqEVX1q1X1jqp6e1W9cDb8s6rql6vqj2f/Hh9dK+xUVUer6m1V9frZ+0dV1Vtnfe+pqrrf6BrhYlX1kKp6TVX9QVWtV9WX62+Zuqr6ttlvhLuq6tVVdUx/yxRV1Sur6oNVdddFw3btY2vby2Zt+Per6vHjKv8rAvkcVNXRJP8+yVcneWySZ1XVY8dWBZd0Lsl3dPdjkzwhybfM2uuLk7yxux+d5I2z9zA1L0yyftH7f5Pkh7r7ryc5k+R5Q6qCS3tpkl/s7sck+aJst1/9LZNVVQ9L8oIkJ7v7cUmOJnlm9LdM06uSPHXHsEv1sV+d5NGz/56f5OULqvGyBPL5+NIkf9Ldf9rdn0jys0meMbgm2FV3393dvzN7/dFs/zh8WLbb7E/OJvvJJP/NkALhEqrq4UmeluQVs/eV5MlJXjObRLtlUqrqwUmemOTHk6S7P9HdH4n+lulbSfKZVbWS5P5J7o7+lgnq7jcn+fCOwZfqY5+R5Kd621uSPKSqPmchhV6GQD4fD0vynovev3c2DCatqh6Z5EuSvDXJanffPRv1gSSro+qCS/jhJN+V5MLs/Wcn+Uh3n5u91/cyNY9KspnkJ2a3Wryiqh4Q/S0T1t3vS/KDSd6d7SB+T5I7o7/l4LhUHzvJzCaQwyFVVQ9M8nNJXtTdf3HxuN7+8wv+BAOTUVVPT/LB7r5zdC1wFVaSPD7Jy7v7S5J8PDsuT9ffMjWz+22fke3/oXRDkgfk0y8JhgPhIPSxAvl8vC/JIy56//DZMJikqvqMbIfxn+nu184Gb9x72c7s3w+Oqg928ZVJ/uuqeme2bwt6crbvzX3I7JLKRN/L9Lw3yXu7+62z96/JdkDX3zJlNyX5s+7e7O5PJnlttvtg/S0HxaX62ElmNoF8Pn47yaNnT5+8X7YffPG6wTXBrmb33f54kvXufslFo16X5Jtmr78pyS8suja4lO7+7u5+eHc/Mtt97K909zcm+dUkXzubTLtlUrr7A0neU1VfOBv0lCTviP6WaXt3kidU1f1nvxnubbf6Ww6KS/Wxr0vynNnT1p+Q5J6LLm0fprbP4nNfVdXfy/b9jUeTvLK7bxlbEeyuqv52kl9Lcjp/dS/u/5Dt+8hvS/K5Sd6V5Ou6e+dDMmC4qnpSku/s7qdX1edn+4z5ZyV5W5Jnd/fWwPLgU1TVF2f7QYT3S/KnSZ6b7RMi+lsmq6r+VZKvz/ZfZnlbkm/O9r22+lsmpapeneRJSR6aZCPJ9yb5+ezSx87+B9OPZPsWjL9M8tzuvmNA2Z9CIAcAAIABXLIOAAAAAwjkAAAAMIBADgAAAAMI5AAAADCAQA4AAAADCOQAsCSq6mNXGP/OqnroIpYFAFyZQA4AAAADCOQAsGSq6nOq6s1V9btVdVdVfdUu0/x8Vd1ZVW+vqudfNPxjVXVLVf1eVb2lqlZnwx9VVb9VVaer6n9e5PoAwLISyAFg+XxDkjd09xcn+aIkv7vLNP9td//NJCeTvKCqPns2/AFJ3tLdX5TkzUn+8Wz4S5O8vLvXkty9j7UDwKEhkAPA8vntJM+tqu9LstbdH91lmhdU1e8leUuSRyR59Gz4J5K8fvb6ziSPnL3+yiSvnr3+3/ehZgA4dARyAFgy3f3mJE9M8r4kr6qq51w8vqqelOSmJF8+OxP+tiTHZqM/2d09e30+ycrFs97HsgHg0BHIAWDJVNXnJdno7h9L8ookj98xyYOTnOnuv6yqxyR5wh5m+xtJnjl7/Y1zKxYADjGBHACWz5OS/F5VvS3J12f7/u+L/WKSlapaT/ID2b5s/UpemORbqup0kofNsVYAOLTqr65KAwAAABbFGXIAAAAYQCAHAACAAQRyAAAAGEAgBwAAgAEEcgAAABhAIAcAAIABBHIAAAAYQCAHAACAAf4/y+w4jbY6kiQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1224x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(17, 6))\n",
    "\n",
    "plt.plot(np.arange(0, 100, dtype=int), positions[0:100], 'o', ms=3)\n",
    "plt.title('King Markov visits')\n",
    "plt.xlabel('Island')\n",
    "plt.ylabel('Week')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fbee26f",
   "metadata": {},
   "source": [
    "### R Code 9.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ee0ac142",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/EAAAGDCAYAAABjtJB2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjb0lEQVR4nO3df9hldV0v/PcnJs2fMaZxDT8E9Jks0CdEjlKmmViiklhPGRRiZqLnSHkefc5JK7MTh/KcUo9cKR1UEh8I9MFfNFFKlnLqiMqv+KXGqBCMI2iMgloo8Hn+2Gtyc889Mxucfe9Z3q/Xde1rr/Vd37XWh31d+2Le9/e7vru6OwAAAMDu77sWXQAAAAAwGyEeAAAARkKIBwAAgJEQ4gEAAGAkhHgAAAAYCSEeAAAARkKIB4AFqaqHV9VXq2qPRddyT1TVdVX1tAXe/+qqespO+ozyswWAnRHiAWCOlgbeqjqmqrZU1Y939z919wO7+8453PftVdVVdfSS9jcM7b+8q++5Urr74O7+8E763O2zraoPV9WvrkiBADBHQjwArJCqen6SNyV5Vnd/ZAVu+Y9Jjp+6/5okz03ymXtzseF8AGCBhHgAWAFV9eIkr0vy9O7+30PbAcOo+Jph/8NVdVJV/X1V3VZVH6yqh05d4/iqur6q/rmqXj3DtPY/T/JjVbV22D8yyRVJvjB1zUdW1d8M1/xSVZ1VVXtOHb+uqn6jqq5I8rWlQb6qfqiqPldVxw77L6qqjVV1S1WdV1V7D+2nVtUfLTn3/VX18mU+qx32nf7vrqrHV9XFVXVrVd1UVa9f+tlW1clJnpTkj4cp9n9cE2+oqpuHc6+sqkfv4LMEgN2CEA8A8/fvk/xekiO6++Kd9P3FJC9I8v1J7pPk/0mSqjooyZuT/FKSdUm+N8k+O7nWvyZ5f5Jjhv3jk7xjSZ9K8gdJ9k7yQ0n2S/K7S/ocm+RZSfbs7jv+7cSqQ5N8IMmvdffZVfXU4VrPHWq8Psk5Q/ezk/xCVdVw7tokPzV1fNo96fvGJG/s7gcneWSSdy3t0N2/leR/JTlxmGJ/4nC9Jyf5gUw+y+cm+edlrg8AuxUhHgDm7yeTXJTkyhn6/ml3/2N3/0smgfSQof3nkvx5d/9dd38jye8k6Rmu944kxw+j6z+e5H3TB7t7Y3df0N23d/cXk7x+6DftlO6+YahpqyclOS/J8d29YWj7pSSnd/el3X17klcl+ZGqOiCTEN3DeVv/ez7a3Z9fpuZ70vebSf6Pqnpod3+1uy/a0Yex5LwHJfnBJNXdn+zuzTOeCwALI8QDwPz9+0xGfN+6dXR5B74wtf31JA8ctvdOcsPWA9399cwwctzdf5fkYUl+K8mGJUE8VbVXVZ1TVZuq6tYkZyZ56JLL3JBtvSTJ/16ywNzemYy+b733V4ca9+nuzmQk/djh8C8mOWs7Nc/cN8kLM/lsP1VVn6iqo7bTb+k9/ibJH2eyRsHNVXVaVT14lnMBYJGEeACYv5uSHJHJyPKb7+U1NifZd+tOVd0vyffNeO6ZSV6RbafSJ8nvZzLq/ZhhSvpxmUyxn7bciP9Lkjy8qt4w1fb5JPtP1fiAocZNQ9PZSX6uqvZP8oQk795BzTP17e5ru/vYTB4/+G9Jzh3uu03XZc49pbsfl+SgTP4Q8J92UA8A7BaEeABYAcNU8COSHLkk+M7q3CQ/XVU/WlX3yeS59Z2N6m91SiZT+i9c5tiDknw1yVeqap/MHmRvy2ShvCdX1WuHtrOTvKCqDqmq+2byB4KPdfd1SdLdlyX5UpK3JvlAd395exeftW9VHVdVD+vuu5Js7XPXMl1vSvKIqfP+XVU9oaq+O8nXMlk/YLnzAGC3IsQDwArp7n9K8tRMRpj/4B6ee3WSX8tkmvnmTIL3zUlun+HcW7r7Q8M09aX+S5JDk3wlyV8kec89qOnLmfxx4BlVdVJ3/3WSV2cyar45k4Xmjlly2p8ledrwvjOz9D0yydVV9dVMFrk7ZukjA4M3ZvK5b6mqU5I8OMlbkmzJ5BGAf07yhzPUBAALVcv//xwA2J1V1QMzGXle392fW3A5AMAKMRIPACNRVT9dVfcfnvn+o0xWu79usVUBACtJiAeA8Tg6k8XjPp9kfSZTx02pA4BVxHR6AAAAGAkj8QAAADASQjwAAACMxJpFFzAva9eu7b333nvRZQAAAMA9cs0113ypux+23LHv2BC/9957553vfOeiywAAAIB75DGPecz12ztmOj0AAACMhBAPAAAAIyHEAwAAwEgI8QAAADASQjwAAACMhBAPAAAAIyHEAwAAwEgI8QAAADASQjwAAACMhBAPAAAAIyHEAwAAwEgI8QAAADASQjwAAACMxJpFFwAAAMDsjjrz+kWXMEobjtt/0SXsEkbiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkRDiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkRDiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkZhbiK+q06vq5qq6aqrtnVV1+fC6rqouH9oPqKp/mTr2J1PnPK6qrqyqjVV1SlXVvGoGAACA3dmaOV777Un+OMk7tjZ09y9s3a6q1yX5ylT/z3T3Ictc59QkL0rysSTnJzkyyV/u+nIBAABg9za3kfjuvjDJLcsdG0bTn5vk7B1do6rWJXlwd1/U3Z3JHwSes4tLBQAAgFFY1DPxT0pyU3dfO9V2YFVdVlUfqaonDW37JLlxqs+NQ9uyquqEqrq4qi7esmXLrq8aAAAAFmhRIf7Y3H0UfnOSh3f3Y5O8PMmfVdWD7+lFu/u07j6suw9bu3btLioVAAAAdg/zfCZ+WVW1JsnPJnnc1rbuvj3J7cP2JVX1mSQ/kGRTkn2nTt93aAMAAIBVZxEj8U9L8qnu/rdp8lX1sKraY9h+RJL1ST7b3ZuT3FpVhw/P0R+f5P0LqBkAAAAWbp4/MXd2ko8meVRV3VhVLxwOHZNtF7R7cpIrhp+cOzfJS7p766J4/yHJW5NsTPKZWJkeAACAVWpu0+m7+9jttP/yMm3vTvLu7fS/OMmjd2lxAAAAMEKLWtgOAAAAuIeEeAAAABgJIR4AAABGQogHAACAkRDiAQAAYCTmtjo9AACweh115vWLLmGUNhy3/6JLYDdnJB4AAABGQogHAACAkRDiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkRDiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkRDiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkRDiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkRDiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkRDiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkRDiAQAAYCTmFuKr6vSqurmqrppq+92q2lRVlw+vZ04de1VVbayqT1fV06fajxzaNlbVK+dVLwAAAOzu5jkS//YkRy7T/obuPmR4nZ8kVXVQkmOSHDyc8+aq2qOq9kjypiTPSHJQkmOHvgAAALDqrJnXhbv7wqo6YMbuRyc5p7tvT/K5qtqY5PHDsY3d/dkkqapzhr7X7Op6AQAAYHe3iGfiT6yqK4bp9muHtn2S3DDV58ahbXvty6qqE6rq4qq6eMuWLbu6bgAAAFiolQ7xpyZ5ZJJDkmxO8rpdefHuPq27D+vuw9auXbvzEwAAAGBE5jadfjndfdPW7ap6S5INw+6mJPtNdd13aMsO2gEAAGBVWdGR+KpaN7X7M0m2rlx/XpJjquq+VXVgkvVJPp7kE0nWV9WBVXWfTBa/O28lawYAAIDdxdxG4qvq7CRPSfLQqroxyWuSPKWqDknSSa5L8uIk6e6rq+pdmSxYd0eSl3b3ncN1TkzygSR7JDm9u6+eV80AAACwO5vn6vTHLtP8th30PznJycu0n5/k/F1YGgAAAIzSIlanBwAAAO4FIR4AAABGQogHAACAkRDiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkZjb78QDAMDu7Kgzr190CaO14bj9F10CrFpG4gEAAGAkhHgAAAAYCSEeAAAARkKIBwAAgJEQ4gEAAGAkhHgAAAAYCSEeAAAARkKIBwAAgJEQ4gEAAGAkhHgAAAAYCSEeAAAARkKIBwAAgJEQ4gEAAGAkhHgAAAAYCSEeAAAARkKIBwAAgJEQ4gEAAGAkhHgAAAAYCSEeAAAARkKIBwAAgJEQ4gEAAGAkhHgAAAAYCSEeAAAARkKIBwAAgJEQ4gEAAGAk5hbiq+r0qrq5qq6aavvDqvpUVV1RVe+tqj2H9gOq6l+q6vLh9SdT5zyuqq6sqo1VdUpV1bxqBgAAgN3ZPEfi357kyCVtFyR5dHf/n0n+Mcmrpo59prsPGV4vmWo/NcmLkqwfXkuvCQAAAKvC3EJ8d1+Y5JYlbR/s7juG3YuS7Luja1TVuiQP7u6LuruTvCPJc+ZQLgAAAOz2FvlM/K8k+cup/QOr6rKq+khVPWlo2yfJjVN9bhzallVVJ1TVxVV18ZYtW3Z9xQAAALBACwnxVfVbSe5IctbQtDnJw7v7sUlenuTPqurB9/S63X1adx/W3YetXbt21xUMAAAAu4E1K33DqvrlJEclOWKYIp/uvj3J7cP2JVX1mSQ/kGRT7j7lft+hDQAAAFadFR2Jr6ojk/znJM/u7q9PtT+sqvYYth+RyQJ2n+3uzUlurarDh1Xpj0/y/pWsGQAAAHYXcxuJr6qzkzwlyUOr6sYkr8lkNfr7Jrlg+KW4i4aV6J+c5Peq6ptJ7kryku7euijef8hkpfv7ZfIM/fRz9AAAALBqzC3Ed/exyzS/bTt9353k3ds5dnGSR+/C0gAAAGCUdjqdvqp+vqoeNGz/dlW9p6oOnX9pAAAAwLRZnol/dXffVlU/luRpmYymnzrfsgAAAIClZgnxdw7vz0pyWnf/RZL7zK8kAAAAYDmzhPhNVfU/k/xCkvOr6r4zngcAAADsQrOE8ecm+UCSp3f3l5M8JMl/mmdRAAAAwLZmCfHHdvd7uvvaJBl+u/2p8y0LAAAAWGqWn5j7v6rqX7v7rCSpqjcl+Z75lgUAAAAsNVOIT3JeVd2V5MgkX+7uF863LAAAAGCp7Yb4qnrI1O6vJnlfkr9P8l+q6iHdfcucawMAAACm7Ggk/pIknaSm3p81vDrJI+ZeHQAAAPBvthviu/vAlSwEAAAA2LGdPhNfVfdP8vIkD+/uE6pqfZJHdfeGuVcHALDKHXXm9YsuYbQ2HLf/oksA2OVm+Ym5P03yjSQ/OuxvSvJf51YRAAAAsKxZQvwju/u/J/lmknT31zN5Ph4AAABYQbOE+G9U1f0yWcwuVfXIJLfPtSoAAABgG7P8TvxrkvxVkv2q6qwkT0zyy/MsCgAAANjWTkN8d19QVZcmOTyTafQv6+4vzb0yAAAA4G52Op2+qirJM5I8bliR/v5V9fi5VwYAAADczSzPxL85yY8kOXbYvy3Jm+ZWEQAAALCsWZ6Jf0J3H1pVlyVJd2+pqvvMuS4AAABgiVlG4r9ZVXvkW6vTPyzJXXOtCgAAANjGLCH+lCTvTfL9VXVykr9L8vtzrQoAAADYxiyr059VVZckOSKT1emf092fnHtlAAAAwN3sNMRX1UlJLkzy9u7+2vxLAgAAAJYzy3T6z2ayMv3FVfXxqnpdVR0957oAAACAJXYa4rv7T7v7V5L8RJIzk/z88A4AAACsoFmm0781yUFJbkryv5L8XJJL51wXAAAAsMQs0+m/L8keSb6c5JYkX+ruO+ZZFAAAALCtWVan/5kkqaofSvL0JH9bVXt0977zLg4AAAD4llmm0x+V5ElJnpxkzyR/k8m0egAAAGAF7TTEJzkyk9D+xu7+/JzrAQAAALZjlun0J65EIQAAAMCOzbKwHQAAALAbmGuIr6rTq+rmqrpqqu0hVXVBVV07vK8d2quqTqmqjVV1RVUdOnXO84f+11bV8+dZMwAAAOyuthviq+pDw/t/+zau//ZMnqmf9sokH+ru9Uk+NOwnyTOSrB9eJyQ5dbj/Q5K8JskTkjw+yWu2Bn8AAABYTXY0Er+uqn40ybOr6rFVdej0a5aLd/eFmfy2/LSjk5wxbJ+R5DlT7e/oiYuS7FlV6zL5WbsLuvuW7t6S5IJs+4cBAAAA+I63o4XtfifJq5Psm+T1S451kqfey3vu1d2bh+0vJNlr2N4nyQ1T/W4c2rbXvo2qOiGTUfysW7fuXpYHAAAAu6fthvjuPjfJuVX16u4+aR437+6uqt6F1zstyWlJcvDBB++y6wIAAMDuYJafmDupqp6d5MlD04e7e8O3cc+bqmpdd28epsvfPLRvSrLfVL99h7ZNSZ6ypP3D38b9AQAAYJR2ujp9Vf1BkpcluWZ4vayqfv/buOd5SbauMP/8JO+faj9+WKX+8CRfGabdfyDJT1XV2mFBu58a2gAAAGBV2elIfJJnJTmku+9Kkqo6I8llSX5zZydW1dmZjKI/tKpuzGSV+dcmeVdVvTDJ9UmeO3Q/P8kzk2xM8vUkL0iS7r6lqk5K8omh3+9199LF8gAAAOA73iwhPkn2zLdWmf/eWS/e3cdu59ARy/TtJC/dznVOT3L6rPcFAFbWUWdev+gSRmvDcfsvugQARmSWEP8HSS6rqr9NUpk8G//KHZ8CAAAA7GqzLGx3dlV9OMm/G5p+o7u/MNeqAAAAgG3MNJ1+WGDuvDnXAgAAAOzATlenBwAAAHYPQjwAAACMxA5DfFXtUVWfWqliAAAAgO3bYYjv7juTfLqqHr5C9QAAAADbMcvCdmuTXF1VH0/yta2N3f3suVUFAAAAbGOWEP/quVcBAAAA7NQsvxP/karaP8n67v7rqrp/kj3mXxoAAAAwbaer01fVi5Kcm+R/Dk37JHnfHGsCAAAAljHLT8y9NMkTk9yaJN19bZLvn2dRAAAAwLZmeSb+9u7+RlUlSapqTZKea1UAMAdHnXn9oksYrQ3H7b/oEgCAzDYS/5Gq+s0k96uqn0zy/yX58/mWBQAAACw1S4h/ZZIvJrkyyYuTnJ/kt+dZFAAAALCtWVanv6uqzkjysUym0X+6u02nBwAAgBW20xBfVc9K8idJPpOkkhxYVS/u7r+cd3EAAADAt8yysN3rkvxEd29Mkqp6ZJK/SCLEAwAAwAqa5Zn427YG+MFnk9w2p3oAAACA7djuSHxV/eyweXFVnZ/kXZk8E//zST6xArUBAAAAU3Y0nf6np7ZvSvLjw/YXk9xvbhUBAAAAy9puiO/uF6xkIQAAAMCOzbI6/YFJfi3JAdP9u/vZ8ysLAAAAWGqW1enfl+RtSf48yV1zrQYAAADYrllC/L929ylzrwQAAADYoVlC/Bur6jVJPpjk9q2N3X3p3KoCAAAAtjFLiH9MkucleWq+NZ2+h30AAABghcwS4n8+ySO6+xvzLgYAAADYvu+aoc9VSfaccx0AAADATswyEr9nkk9V1Sdy92fi/cQcAAAArKBZQvxr5l4FAAAAsFM7DfHd/ZGVKAQAAADYsZ2G+Kq6LZPV6JPkPkm+O8nXuvvB8ywMAAAAuLtZRuIftHW7qirJ0UkOv7c3rKpHJXnnVNMjkvxOJs/evyjJF4f23+zu84dzXpXkhUnuTPLr3f2Be3t/AAAAGKtZVqf/Nz3xviRPv7c37O5Pd/ch3X1Ikscl+XqS9w6H37D12FSAPyjJMUkOTnJkkjdX1R739v4AAAAwVrNMp//Zqd3vSnJYkn/dRfc/Islnuvv6ySD/so5Ock53357kc1W1Mcnjk3x0F9UAAAAAozDL6vQ/PbV9R5LrMgnWu8IxSc6e2j+xqo5PcnGSV3T3liT7JLloqs+NQxvAd6Sjzrx+0SWM0obj9l90CQAAczfLM/EvmMeNq+o+SZ6d5FVD06lJTspkEb2Tkrwuya/cw2uekOSEJFm3bt0uqxUAAAB2B9sN8VX1Ozs4r7v7pG/z3s9Icml33zRc8Kape78lyYZhd1OS/abO23doW66o05KcliQHH3xwL9cHAAAAxmpHC9t9bZlXMlkl/jd2wb2PzdRU+qqaHjr/mSRXDdvnJTmmqu5bVQcmWZ/k47vg/gAAADAq2x2J7+7Xbd2uqgcleVmSFyQ5J5Op7vdaVT0gyU8mefFU83+vqkMymU5/3dZj3X11Vb0ryTWZPJP/0u6+89u5PwAAAIzRDp+Jr6qHJHl5kl9KckaSQ4fF5r4t3f21JN+3pO15O+h/cpKTv937AgAAwJjt6Jn4P0zys5k8Y/6Y7v7qilUFAAAAbGNHz8S/IsneSX47yeer6tbhdVtV3boy5QEAAABb7eiZ+B0FfAAAAGCFCeoAAAAwEkI8AAAAjIQQDwAAACMhxAMAAMBICPEAAAAwEkI8AAAAjIQQDwAAACMhxAMAAMBICPEAAAAwEmsWXQAwTkedef2iSxilDcftv+gSAAAYMSPxAAAAMBJCPAAAAIyEEA8AAAAjIcQDAADASAjxAAAAMBJCPAAAAIyEEA8AAAAjIcQDAADASAjxAAAAMBJCPAAAAIyEEA8AAAAjIcQDAADASAjxAAAAMBJCPAAAAIyEEA8AAAAjIcQDAADASAjxAAAAMBJCPAAAAIyEEA8AAAAjIcQDAADASCwsxFfVdVV1ZVVdXlUXD20PqaoLqura4X3t0F5VdUpVbayqK6rq0EXVDQAAAIuy6JH4n+juQ7r7sGH/lUk+1N3rk3xo2E+SZyRZP7xOSHLqilcKAAAAC7boEL/U0UnOGLbPSPKcqfZ39MRFSfasqnULqA8AAAAWZpEhvpN8sKouqaoThra9unvzsP2FJHsN2/skuWHq3BuHNgAAAFg11izw3j/W3Zuq6vuTXFBVn5o+2N1dVX1PLjj8MeCEJFm3zkA9AAAA31kWNhLf3ZuG95uTvDfJ45PctHWa/PB+89B9U5L9pk7fd2hbes3Tuvuw7j5s7dq18ywfAAAAVtxCQnxVPaCqHrR1O8lPJbkqyXlJnj90e36S9w/b5yU5flil/vAkX5madg8AAACrwqKm0++V5L1VtbWGP+vuv6qqTyR5V1W9MMn1SZ479D8/yTOTbEzy9SQvWPmSAQAAYLEWEuK7+7NJfniZ9n9OcsQy7Z3kpStQGgAAAOy2drefmAMAAAC2Q4gHAACAkRDiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkRDiAQAAYCQW8jvxsFKOOvP6RZcwShuO23/RJQAAAMswEg8AAAAjIcQDAADASAjxAAAAMBJCPAAAAIyEEA8AAAAjIcQDAADASAjxAAAAMBJCPAAAAIyEEA8AAAAjIcQDAADASAjxAAAAMBJCPAAAAIyEEA8AAAAjIcQDAADASAjxAAAAMBJCPAAAAIyEEA8AAAAjIcQDAADASAjxAAAAMBJCPAAAAIyEEA8AAAAjIcQDAADASAjxAAAAMBJCPAAAAIyEEA8AAAAjseIhvqr2q6q/raprqurqqnrZ0P67VbWpqi4fXs+cOudVVbWxqj5dVU9f6ZoBAABgd7BmAfe8I8kruvvSqnpQkkuq6oLh2Bu6+4+mO1fVQUmOSXJwkr2T/HVV/UB337miVQMAAMCCrfhIfHdv7u5Lh+3bknwyyT47OOXoJOd09+3d/bkkG5M8fv6VAgAAwO5loc/EV9UBSR6b5GND04lVdUVVnV5Va4e2fZLcMHXajdlx6AcAAIDvSAsL8VX1wCTvTvIfu/vWJKcmeWSSQ5JsTvK6e3HNE6rq4qq6eMuWLbuyXAAAAFi4hYT4qvruTAL8Wd39niTp7pu6+87uvivJW/KtKfObkuw3dfq+Q9s2uvu07j6suw9bu3btcl0AAABgtBaxOn0leVuST3b366fa1011+5kkVw3b5yU5pqruW1UHJlmf5OMrVS8AAADsLhaxOv0TkzwvyZVVdfnQ9ptJjq2qQ5J0kuuSvDhJuvvqqnpXkmsyWdn+pVamBwAAYDVa8RDf3X+XpJY5dP4Ozjk5yclzKwoAAABGYKGr0wMAAACzE+IBAABgJIR4AAAAGIlFLGzHEkedef2iSxilDcftv+gSAAAAVpSReAAAABgJIR4AAABGQogHAACAkRDiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkRDiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkRDiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkRDiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkRDiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkRDiAQAAYCSEeAAAABgJIR4AAABGQogHAACAkRDiAQAAYCSEeAAAABiJ0YT4qjqyqj5dVRur6pWLrgcAAABW2ihCfFXtkeRNSZ6R5KAkx1bVQYutCgAAAFbWKEJ8kscn2djdn+3ubyQ5J8nRC64JAAAAVtRYQvw+SW6Y2r9xaAMAAIBVo7p70TXsVFX9XJIju/tXh/3nJXlCd5+4pN8JSU4Ydh+V5NMrWihj8dAkX1p0EbBgvgesdr4D4HsAu/N3YP/ufthyB9asdCX30qYk+03t7zu03U13n5bktJUqinGqqou7+7BF1wGL5HvAauc7AL4HMNbvwFim038iyfqqOrCq7pPkmCTnLbgmAAAAWFGjGInv7juq6sQkH0iyR5LTu/vqBZcFAAAAK2oUIT5Juvv8JOcvug6+I3jkAnwPwHcAfA9glN+BUSxsBwAAAIznmXgAAABY9YR4Vo2q2q+q/raqrqmqq6vqZYuuCRahqvaoqsuqasOia4FFqKo9q+rcqvpUVX2yqn5k0TXBSqqq/3v4t9BVVXV2VX3PomuCeauq06vq5qq6aqrtIVV1QVVdO7yvXWSNsxLiWU3uSPKK7j4oyeFJXlpVBy24JliElyX55KKLgAV6Y5K/6u4fTPLD8X1gFamqfZL8epLDuvvRmSwafcxiq4IV8fYkRy5pe2WSD3X3+iQfGvZ3e0I8q0Z3b+7uS4ft2zL5R9s+i60KVlZV7ZvkWUneuuhaYBGq6nuTPDnJ25Kku7/R3V9eaFGw8tYkuV9VrUly/ySfX3A9MHfdfWGSW5Y0H53kjGH7jCTPWcma7i0hnlWpqg5I8tgkH1twKbDS/keS/5zkrgXXAYtyYJIvJvnT4bGSt1bVAxZdFKyU7t6U5I+S/FOSzUm+0t0fXGxVsDB7dffmYfsLSfZaZDGzEuJZdarqgUneneQ/dveti64HVkpVHZXk5u6+ZNG1wAKtSXJoklO7+7FJvpaRTJ+EXWF45vfoTP6gtXeSB1TVcYutChavJz/bNoqfbhPiWVWq6rszCfBndfd7Fl0PrLAnJnl2VV2X5JwkT62qMxdbEqy4G5Pc2N1bZ2Kdm0moh9XiaUk+191f7O5vJnlPkh9dcE2wKDdV1bokGd5vXnA9MxHiWTWqqjJ5BvKT3f36RdcDK627X9Xd+3b3AZksYvQ33W30hVWlu7+Q5IaqetTQdESSaxZYEqy0f0pyeFXdf/i30RGxuCOr13lJnj9sPz/J+xdYy8yEeFaTJyZ5Xiajj5cPr2cuuigAVtyvJTmrqq5IckiS319sObByhlko5ya5NMmVmeSB0xZaFKyAqjo7yUeTPKqqbqyqFyZ5bZKfrKprM5ml8tpF1jirmkz9BwAAAHZ3RuIBAABgJIR4AAAAGAkhHgAAAEZCiAcAAICREOIBAABgJIR4AFjFquqrOzl+XVU9dCXuBQDsnBAPAAAAIyHEAwCpqnVVdWFVXV5VV1XVk5bp876quqSqrq6qE6bav1pVJ1fVP1TVRVW119B+YFV9tKqurKr/upL/PQDwnUqIBwCS5BeTfKC7D0nyw0kuX6bPr3T345IcluTXq+r7hvYHJLmou384yYVJXjS0vzHJqd39mCSb51g7AKwaQjwAkCSfSPKCqvrdJI/p7tuW6fPrVfUPSS5Ksl+S9UP7N5JsGLYvSXLAsP3EJGcP2//vHGoGgFVHiAcA0t0XJnlykk1J3l5Vx08fr6qnJHlakh8ZRtwvS/I9w+FvdncP23cmWTN96TmWDQCrjhAPAKSq9k9yU3e/Jclbkxy6pMv3JtnS3V+vqh9McvgMl/37JMcM27+0y4oFgFVMiAcAkuQpSf6hqi5L8guZPM8+7a+SrKmqTyZ5bSZT6nfmZUleWlVXJtlnF9YKAKtWfWv2GwAAALA7MxIPAAAAIyHEAwAAwEgI8QAAADASQjwAAACMhBAPAAAAIyHEAwAAwEgI8QAAADASQjwAAACMxP8PRa85ZEVzIgkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1224x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(17, 6))\n",
    "\n",
    "plt.hist(positions, rwidth=0.9)\n",
    "plt.title('King Markov visits')\n",
    "plt.xlabel('Island')\n",
    "plt.ylabel('Number of weeks')\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c271859b",
   "metadata": {},
   "source": [
    "### R Code 9.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c3bc60f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def multivariate_normal(D, T=1000):\n",
    "    \"\"\" Return the Gaussian Multivariate standart D-dimensional \n",
    "    \n",
    "    Parameter:\n",
    "    ===========\n",
    "    D: integer\n",
    "        D is the dimension of the Gaussian Multivariate\n",
    "        \n",
    "    T: Integer; default=1000\n",
    "        Qty samples\n",
    "    \"\"\"\n",
    "    means = np.repeat(0, D) \n",
    "    diag = np.diag(np.repeat(1, D))\n",
    "    return np.random.multivariate_normal(means, diag, T)\n",
    "\n",
    "\n",
    "def rad_dist(Y):\n",
    "    \"\"\" Return the distance radial from the mode\n",
    "    \n",
    "    Parameter:\n",
    "    ===========\n",
    "    Y: Multivariate Normal Starndart\n",
    "    \"\"\"\n",
    "    return np.sqrt(np.sum(np.power(Y, 2)))\n",
    "\n",
    "\n",
    "def get_radial_distances(D, T=1000):\n",
    "    \"\"\" Return the array of distances of Gaussian Multivariate standart D-dimensional from mode \n",
    "    \n",
    "    Parameter:\n",
    "    ===========\n",
    "    D: integer\n",
    "        D is the dimension of the Gaussian Multivariate\n",
    "        \n",
    "    T: Integer; default=1000\n",
    "        Qty samples\n",
    "    \"\"\"\n",
    "    Y = multivariate_normal(D, T)\n",
    "    return [rad_dist(Y[i]) for i in range(T)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e860e730",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9kAAAFzCAYAAAAqrMpLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfHUlEQVR4nO3df5RnZ10f8PeHXaMEUnYV9CTZ1ERcf6RAgzsm/ZFasOAJ9DSRCjap9kBrjbRGbbGexupJMW2t0ta2nqaUqFRsiyFSoGu7Gqng0ViIOwvhR5JG1oDNJlYq7NakUSDy6R9zg1+H2ZnvJs/8fr3OmTP3Pvf53vuZud9zd9/z3Pt8q7sDAAAAPHFP2uwCAAAAYKcQsgEAAGAQIRsAAAAGEbIBAABgECEbAAAABhGyAQAAYJC9m3Xg/fv393nnnbdZhwcAAIDH7e677/6d7n7G8vZNC9nnnXde3vSmN23W4QEAAOBxe/azn/2bK7W7XRwAAAAGEbIBAABgECEbAAAABhGyAQAAYBAhGwAAAAYRsgEAAGAQIRsAAAAGEbIBAABgECEbAAAABhGyAQAAYBAhGwAAAAYRsgEAAGCQuUJ2VV1RVfdW1fGqun6F7X+8qt5ZVe+tqvdX1YvHlwoAAABb2961OlTVniQ3JXlhkhNJjlbV4e6+e6bb9ye5tbtfW1UXJzmS5MJ1qHdzLC6u3L6wsLF1AAAAsKXNM5J9aZLj3X1fd38yyS1JrlrWp5P8sWn5aUkeHFciAAAAbA9rjmQnOT/J/TPrJ5JctqzPq5P8QlV9R5KnJHnBSjuqqmuTXJsk55577pnWCgAAAFvaqInPrknyk919IMmLk/yHqvqsfXf3zd290N0L+/fvH3RoAAAA2BrmCdkPJLlgZv3A1DbrW5LcmiTd/a4kn5fk6SMKBAAAgO1inpB9NMnBqrqoqs5KcnWSw8v6/K8kfyFJquorsxSy/8/IQgEAAGCrWzNkd/ejSa5LcluSe7I0i/hdVXVjVV05dfvuJN9aVe9L8tNJXtHdvV5FAwAAwFY0z8Rn6e4jWfpYrtm2G2aW707yZ8eWBgAAANvLqInPAAAAYNcTsgEAAGAQIRsAAAAGEbIBAABgECEbAAAABhGyAQAAYBAhGwAAAAYRsgEAAGAQIRsAAAAGEbIBAABgECEbAAAABhGyAQAAYBAhGwAAAAYRsgEAAGAQIRsAAAAGEbIBAABgECEbAAAABhGyAQAAYBAhGwAAAAYRsgEAAGAQIRsAAAAGEbIBAABgECEbAAAABhGyAQAAYBAhGwAAAAYRsgEAAGAQIRsAAAAGEbIBAABgECEbAAAABtk7T6equiLJv06yJ8mPd/cPLdv+L5M8f1o9O8kXdve+gXVubYuLn922sLDxdQAAALCp1gzZVbUnyU1JXpjkRJKjVXW4u+9+rE93/92Z/t+R5LnrUCsAAABsafPcLn5pkuPdfV93fzLJLUmuWqX/NUl+ekRxAAAAsJ3ME7LPT3L/zPqJqe2zVNUXJ7koyTueeGkAAACwvYye+OzqJG/u7j9YaWNVXVtVi1W1ePLkycGHBgAAgM01T8h+IMkFM+sHpraVXJ1VbhXv7pu7e6G7F/bv3z9/lQAAALANzBOyjyY5WFUXVdVZWQrSh5d3qqqvSLI/ybvGlggAAADbw5ohu7sfTXJdktuS3JPk1u6+q6purKorZ7peneSW7u71KRUAAAC2trk+J7u7jyQ5sqzthmXrrx5XFgAAAGw/oyc+AwAAgF1LyAYAAIBBhGwAAAAYRMgGAACAQYRsAAAAGETIBgAAgEGEbAAAABhEyAYAAIBBhGwAAAAYRMgGAACAQYRsAAAAGGTvZhew4y0urty+sLCxdQAAALDujGQDAADAIEI2AAAADCJkAwAAwCBCNgAAAAwiZAMAAMAgQjYAAAAMImQDAADAIEI2AAAADCJkAwAAwCBCNgAAAAwiZAMAAMAgQjYAAAAMImQDAADAIEI2AAAADCJkAwAAwCBCNgAAAAwiZAMAAMAgQjYAAAAMMlfIrqorqureqjpeVdefps83VtXdVXVXVb1xbJkAAACw9e1dq0NV7UlyU5IXJjmR5GhVHe7uu2f6HEzyvUn+bHefrKovXK+CAQAAYKuaZyT70iTHu/u+7v5kkluSXLWsz7cmuam7TyZJd390bJkAAACw9c0Tss9Pcv/M+ompbdaXJfmyqvrVqnp3VV2x0o6q6tqqWqyqxZMnTz6+igEAAGCLGjXx2d4kB5M8L8k1SX6sqvYt79TdN3f3Qncv7N+/f9ChAQAAYGuYJ2Q/kOSCmfUDU9usE0kOd/enuvvDSX49S6EbAAAAdo15QvbRJAer6qKqOivJ1UkOL+vztiyNYqeqnp6l28fvG1cmAAAAbH1rhuzufjTJdUluS3JPklu7+66qurGqrpy63ZbkY1V1d5J3Jvme7v7YehUNAAAAW9GaH+GVJN19JMmRZW03zCx3kldNXwAAALArjZr4DAAAAHY9IRsAAAAGEbIBAABgECEbAAAABplr4jNYy7FTx9bsc2jfoQ2oBAAAYPMYyQYAAIBBhGwAAAAYRMgGAACAQYRsAAAAGETIBgAAgEHMLg4AALCF3HHikVW3X3bg7A2qhMfDSDYAAAAMImQDAADAIEI2AAAADCJkAwAAwCBCNgAAAAwiZAMAAMAgQjYAAAAMImQDAADAIEI2AAAADCJkAwAAwCBCNgAAAAyyd7MLANhI5zx4+6rbHzrv8g2qBACAnchINgAAAAxiJBtgFUa+AQA4E0I2AADANnXHiUdW3X7ZgbM3qBIe43ZxAAAAGETIBgAAgEGEbAAAABhkrpBdVVdU1b1Vdbyqrl9h+yuq6v9U1Z3T198cXyoAAABsbWtOfFZVe5LclOSFSU4kOVpVh7v77mVd39Td161DjQAAALAtzDOSfWmS4919X3d/MsktSa5a37IAAABg+5knZJ+f5P6Z9RNT23LfUFXvr6o3V9UFQ6oDAACAbWTUxGc/m+TC7n5OkrcnecNKnarq2qparKrFkydPDjo0AAAAbA3zhOwHksyOTB+Y2j6juz/W3Z+YVn88yaGVdtTdN3f3Qncv7N+///HUCwAAAFvWPCH7aJKDVXVRVZ2V5Ookh2c7VNW5M6tXJrlnXIkAAACwPaw5u3h3P1pV1yW5LcmeJK/v7ruq6sYki919OMl3VtWVSR5N8vEkr1jHmgEAAGBLWjNkJ0l3H0lyZFnbDTPL35vke8eWBgAAANvLXCGb3efYqWOrbj+0b8XH7gEAAHa1UbOLAwAAwK4nZAMAAMAgQjYAAAAMImQDAADAIEI2AAAADCJkAwAAwCBCNgAAAAwiZAMAAMAgQjYAAAAMImQDAADAIEI2AAAADCJkAwAAwCBCNgAAAAwiZAMAAMAgQjYAAAAMImQDAADAIEI2AAAADCJkAwAAwCB7N7uAXW9x8bPbFhY2vg4AAACeMCPZAAAAMIiRbDbFsVPH1uxzaN+hDagEAABgHCPZAAAAMIiQDQAAAIMI2QAAADCIkA0AAACDCNkAAAAwiJANAAAAgwjZAAAAMIiQDQAAAIPMFbKr6oqqureqjlfV9av0+4aq6qpaGFciAAAAbA9rhuyq2pPkpiQvSnJxkmuq6uIV+p2T5LuS3DG6SAAAANgO5hnJvjTJ8e6+r7s/meSWJFet0O8fJfnhJL8/sD4AAADYNuYJ2ecnuX9m/cTU9hlV9VVJLuju/7bajqrq2qparKrFkydPnnGxAAAAsJU94YnPqupJSX4kyXev1be7b+7uhe5e2L9//xM9NAAAAGwp84TsB5JcMLN+YGp7zDlJnpXkl6rqI0n+VJLDJj8DAABgt5knZB9NcrCqLqqqs5JcneTwYxu7+/9299O7+8LuvjDJu5Nc2d2L61IxAAAAbFFrhuzufjTJdUluS3JPklu7+66qurGqrlzvAgEAAGC72DtPp+4+kuTIsrYbTtP3eU+8LAAAANh+nvDEZwAAAMASIRsAAAAGEbIBAABgECEbAAAABhGyAQAAYBAhGwAAAAYRsgEAAGAQIRsAAAAGEbIBAABgECEbAAAABtm72QUA7CTnPHj7qtsfOu/yDaoEAIDNYCQbAAAABjGSDQAAsM7uOPHIqtsvO3D2BlXCejOSDQAAAIMI2QAAADCIkA0AAACDCNkAAAAwiJANAAAAgwjZAAAAMIiQDQAAAIMI2QAAADCIkA0AAACD7N3sAgCeiHMevH3V7Q+dd/kGVQIAsLXdceKRVbdfduDsDapkZzOSDQAAAIMI2QAAADCIkA0AAACDCNkAAAAwiJANAAAAg5hdfKtbXPzstoWFja8DAACANRnJBgAAgEHmCtlVdUVV3VtVx6vq+hW2v7KqPlBVd1bV7VV18fhSAQAAYGtbM2RX1Z4kNyV5UZKLk1yzQoh+Y3c/u7svSfKaJD8yulAAAADY6uYZyb40yfHuvq+7P5nkliRXzXbo7t+dWX1Kkh5XIgAAAGwP80x8dn6S+2fWTyS5bHmnqvr2JK9KclaSr11pR1V1bZJrk+Tcc88901oBAABgSxs28Vl339Tdz0zy95N8/2n63NzdC929sH///lGHBgAAgC1hnpD9QJILZtYPTG2nc0uSr38CNQEAAMC2NE/IPprkYFVdVFVnJbk6yeHZDlV1cGb1Lyb50LgSAQAAYHtY85ns7n60qq5LcluSPUle3913VdWNSRa7+3CS66rqBUk+leRkkpevZ9EAAACwFc0z8Vm6+0iSI8vabphZ/q7BdQEAAMC2M2ziMwAAANjthGwAAAAYRMgGAACAQYRsAAAAGETIBgAAgEGEbAAAABhEyAYAAIBBhGwAAAAYRMgGAACAQYRsAAAAGETIBgAAgEH2bnYBAAAAbC13nHhk1e2XHTh7gyrZfoxkAwAAwCBCNgAAAAwiZAMAAMAgQjYAAAAMImQDAADAIEI2AAAADCJkAwAAwCBCNgAAAAwiZAMAAMAgQjYAAAAMImQDAADAIEI2AAAADCJkAwAAwCB7N7sAnoDFxZXbFxY2tg4AAACSGMkGAACAYYxk7wLHTh1bs8+hfYc2oBIAAICdzUg2AAAADDJXyK6qK6rq3qo6XlXXr7D9VVV1d1W9v6p+saq+eHypAAAAsLWtGbKrak+Sm5K8KMnFSa6pqouXdXtvkoXufk6SNyd5zehCAQAAYKubZyT70iTHu/u+7v5kkluSXDXbobvf2d2PTKvvTnJgbJkAAACw9c0Tss9Pcv/M+omp7XS+JcnPPZGiAAAAYDsaOrt4VX1zkoUkf/40269Ncm2SnHvuuSMPDQAAAJtunpHsB5JcMLN+YGr7I6rqBUm+L8mV3f2JlXbU3Td390J3L+zfv//x1AsAAABb1jwh+2iSg1V1UVWdleTqJIdnO1TVc5O8LksB+6PjywQAAICtb82Q3d2PJrkuyW1J7klya3ffVVU3VtWVU7d/luSpSX6mqu6sqsOn2R0AAADsWHM9k93dR5IcWdZ2w8zyCwbXBQAAANvOPLeLAwAAAHMQsgEAAGAQIRsAAAAGGfo52QAAALvFHSceWXX7ZQfO3qBK2EqEbLaNY6eOrbr90L5DG1QJAADAyoRsgE1yzoO3r9nnofMu34BKAAAYxTPZAAAAMIiQDQAAAIMI2QAAADCIkA0AAACDCNkAAAAwiJANAAAAgwjZAAAAMIiQDQAAAIMI2QAAADCIkA0AAACDCNkAAAAwiJANAAAAgwjZAAAAMIiQDQAAAIPs3ewCAGad8+Dta/Z56LzLN6ASAAA4c0ayAQAAYBAhGwAAAAYRsgEAAGAQIRsAAAAGEbIBAABgECEbAAAABhGyAQAAYBCfk72TLS4mSfY9fPwzTacuObhZ1QAAAOx4RrIBAABgkLlCdlVdUVX3VtXxqrp+he1fU1XvqapHq+ql48sEAACArW/NkF1Ve5LclORFSS5Ock1VXbys2/9K8ookbxxdIAAAAGwX8zyTfWmS4919X5JU1S1Jrkpy92Mduvsj07ZPr0ONAAAAsC3Mc7v4+Unun1k/MbWdsaq6tqoWq2rx5MmTj2cXAAAAsGVt6MRn3X1zdy9098L+/fs38tAAAACw7uYJ2Q8kuWBm/cDUBgAAAMyYJ2QfTXKwqi6qqrOSXJ3k8PqWBQAAANvPmiG7ux9Ncl2S25Lck+TW7r6rqm6sqiuTpKq+uqpOJHlZktdV1V3rWTQAAABsRfPMLp7uPpLkyLK2G2aWj2bpNnIAAADYtTZ04jMAAADYyeYayYbt6NipY2v2ObTv0AZUAgAA7BZGsgEAAGAQI9nb0GMjtPsePr7i9lOnyggtAADAJjCSDQAAAIMI2QAAADCIkA0AAACDCNkAAAAwiJANAAAAg5hdHGCbOOfB21fd/tB5l29QJQAAnI6RbAAAABjESPYutu/OD/3hylN76fvCwuYUAwAAsAMYyQYAAIBBjGQDAADwuN1x4pFVt1924OwNqmRrELIBAAAmAiNPlNvFAQAAYBAj2cC6WutjpxIfPQUAwM5hJBsAAAAGEbIBAABgECEbAAAABhGyAQAAYBATn/H4LC7+kdV9Dx/PqUsOblIxAAAAW4OQDQAAwIbYDZ9DLmTD5NipY2v2ObTv0AZUAgAAbFdCNivad+eHPqvN7eAAAACrM/EZAAAADGIkG2CHOufB21fd/tB5l29QJQAAu4eRbAAAABjESDYAALCj7YYZrdk65hrJrqorqureqjpeVdevsP1zq+pN0/Y7qurC4ZUCAADAFrfmSHZV7UlyU5IXJjmR5GhVHe7uu2e6fUuSk939pVV1dZIfTvJX1qPgneg3Hj7+R9ZPnarPLO97+Hie+dQv3eiSmNNaH/vlI78AAGB3med28UuTHO/u+5Kkqm5JclWS2ZB9VZJXT8tvTvJvqqq6uwfWCmwRJtQCAICVzXO7+PlJ7p9ZPzG1rdinux9N8n+TfMGIAgEAAGC72NCJz6rq2iTXTqsPP/vZz753I4//OD09ye9sdhFsGucf74Hdzfnf3Zx/vAd2N+eftd4DX7xS4zwh+4EkF8ysH5jaVupzoqr2Jnlako8t31F335zk5jmOuWVU1WJ3L2x2HWwO5x/vgd3N+d/dnH+8B3Y355/H+x6Y53bxo0kOVtVFVXVWkquTHF7W53CSl0/LL03yDs9jAwAAsNusOZLd3Y9W1XVJbkuyJ8nru/uuqroxyWJ3H07yE0n+Q1UdT/LxLAVxAAAA2FXmeia7u48kObKs7YaZ5d9P8rKxpW0Z2+r2doZz/vEe2N2c/93N+cd7YHdz/nlc74FyVzcAAACMMc8z2QAAAMAchOxVVNUVVXVvVR2vqus3ux42VlV9pKo+UFV3VtXiZtfD+quq11fVR6vqgzNtn19Vb6+qD03f929mjayf05z/V1fVA9N14M6qevFm1sj6qaoLquqdVXV3Vd1VVd81tbsG7AKrnH/XgF2iqj6vqn6tqt43vQd+YGq/qKrumPLAm6aJoNlhVjn/P1lVH565Blwy1/7cLr6yqtqT5NeTvDDJiSzNsn5Nd9+9qYWxYarqI0kWutvnI+4SVfU1SR5O8lPd/ayp7TVJPt7dPzT9sW1/d//9zayT9XGa8//qJA939z/fzNpYf1V1bpJzu/s9VXVOkmNJvj7JK+IasOOtcv6/Ma4Bu0JVVZKndPfDVfU5SW5P8l1JXpXkLd19S1X9uyTv6+7XbmatjLfK+X9lkv/a3W8+k/0ZyT69S5Mc7+77uvuTSW5JctUm1wSso+7+5Sx9QsKsq5K8YVp+Q5b+08UOdJrzzy7R3b/V3e+Zlh9Kck+S8+MasCuscv7ZJXrJw9Pq50xfneRrkzwWsFwDdqhVzv/jImSf3vlJ7p9ZPxEX292mk/xCVR2rqms3uxg2zRd1929Ny/87yRdtZjFsiuuq6v3T7eRuFd4FqurCJM9NckdcA3adZec/cQ3YNapqT1XdmeSjSd6e5DeSnOruR6cu8sAOtvz8d/dj14B/Ml0D/mVVfe48+xKy4fQu7+6vSvKiJN8+3UrKLtZLz9d4xmZ3eW2SZya5JMlvJfkXm1oN666qnprkPyf5O939u7PbXAN2vhXOv2vALtLdf9DdlyQ5kKW7Wr9icytiIy0//1X1rCTfm6X3wVcn+fwkcz0uJGSf3gNJLphZPzC1sUt09wPT948meWuWLrbsPr89Pav32DN7H93kethA3f3b0z+6n07yY3Ed2NGm5/D+c5L/1N1vmZpdA3aJlc6/a8Du1N2nkrwzyZ9Osq+q9k6b5IFdYOb8XzE9StLd/Ykk/z5zXgOE7NM7muTgNKPgWUmuTnJ4k2tig1TVU6aJT1JVT0nydUk+uPqr2KEOJ3n5tPzyJP9lE2thgz0WriYvievAjjVNevMTSe7p7h+Z2eQasAuc7vy7BuweVfWMqto3LT85S5Mf35OlsPXSqZtrwA51mvP/P2f+yFpZeh5/rmuA2cVXMX1Mw79KsifJ67v7n2xuRWyUqvqSLI1eJ8neJG90/ne+qvrpJM9L8vQkv53kHyZ5W5Jbk/zxJL+Z5Bu72+RYO9Bpzv/zsnSbaCf5SJJvm3k+lx2kqi5P8itJPpDk01PzP8jSc7muATvcKuf/mrgG7ApV9ZwsTWy2J0sDkbd2943T/wlvydKtwu9N8s3TqCY7yCrn/x1JnpGkktyZ5JUzE6Sdfn9CNgAAAIzhdnEAAAAYRMgGAACAQYRsAAAAGETIBgAAgEGEbAAAABhEyAZgR6mqP6iqO6vqg1X1s4997uUZvP6XqmphWj6y1uur6iNV9fTTtH9g+rq7qv5xVX3etO28qnrzKvvcV1V/+0zqHqmqXlZV91TVO9f5OD9ZVY9U1Tkzbf+qqnql3+kq+3l1Vf299akSAM6MkA3ATvN73X1Jdz8ryceTfPvj3VF3v7i7Tz2BWp7f3c9OcmmSL0nyumm/D3b3S1d53b4kmxayk3xLkm/t7ufPNlbV3nU41vEkV037f1KSr03ywDocBwA2hJANwE72riTnJ0lVXVpV76qq91bV/6iqL5/an1xVt0wjt29N8uTHXjw7Sl1Vb6uqY1V1V1VdeyZFdPfDSV6Z5Our6vOr6sKq+uC03z9RVb82jb6/v6oOJvmhJM+c2v5ZVT21qn6xqt4zjYw/FkovnOr+samuX6iqJ0/bvrSq/ntVvW963TOn9u+pqqPTsX5gea1VdUOSy5P8xHTsV1TV4ap6R5JfnOp/2/T6d1fVc6bXvbqq3lBVv1JVv1lVf7mqXjPV+/NV9Tmn+fXckuSvTMvPS/KrSR6dqedV010JH6yqvzPT/n1V9etVdXuSL59pf+Z0vGNTLV9xJucKAJ4oIRuAHamq9iT5C0kOT03/M8mf6+7nJrkhyQ9O7X8rySPd/ZVJ/mGSQ6fZ5d/o7kNJFpJ8Z1V9wZnU092/m+TDSQ4u2/TKJP+6uy+Z9n0iyfVJfmMakf+eJL+f5CXd/VVJnp/kX1RVTa8/mOSm7v4TSU4l+Yap/T9N7X8yyZ9J8ltV9XVT/0uTXJLkUFV9zbI6b0yymOSbpmMnyVcleWl3//kkP5Dkvd39nCT/IMlPzbz8mVkaib4yyX9M8s5pJP/3kvzF0/xqfj3JM6pqf5JrshS6kyRVdSjJX09yWZI/leRbq+q5U/vV08/w4iRfPbO/m5N8x3Su/l6Sf3ua4wLAuliP274AYDM9uaruzNII9j1J3j61Py3JG6aR4k7y2Mjq1yT50STp7vdX1ftPs9/vrKqXTMsXZCmsfuwMa6sV2t6V5Puq6kCSt3T3h/4wP/+R1/3gFIg/Pf1sXzRt+3B33zktH0ty4fSM8/nd/dbp5/r9JJlC9tclee/U/6nTz/HLa9T99u7++LR8eaYg393vqKovqKo/Nm37ue7+VFV9IMmeJD8/tX8gyYWr7P8tWQrNlyX5tpn2y5O8tbv/31T/W5L8uSwNEry1ux+Z2g9P35+apT8o/MzM7/Bz1/jZAGAoIRuAneb3uvuSqjo7yW1Zeib7R5P8oyyNrL6kqi5M8kvz7rCqnpfkBUn+dHc/UlW/lOTzzqSoKfhemKWR26c91t7db6yqO7I00nukqr4tyX3LXv5NSZ6R5NAUYj8yc/xPzPT7g8zc7r5SGUn+aXe/7kxqT/L/5uz3iSTp7k9X1ae6u6f2T2f1/3O8KUt/IHjD9NozLO8znpTk1HRXAABsCreLA7AjTaOc35nku6cJu56WP5xQ6xUzXX85yV9Nkqp6VpLnrLC7pyU5OQXsr8jSrctzm0ZY/22St3X3yWXbviTJfd39o0n+y3T8h5KcM9PtaUk+OgXs5yf54tWO190PJTlRVV8/HeNzZ/7o8DemelJV51fVF57Jz5LkV7IU+h/748PvTLfCP27d/ZtJvi+ffWv3r2TpOfazq+opSV4ytf3y1P7k6Y8Xf2naz+8m+XBVvWyqr6rqTz6R2gDgTBnJBmDH6u73Trd/X5PkNVm6Xfz7k/y3mW6vTfLvq+qeLN1efmyFXf18kldOfe5N8u45S3jn9Oz0k5K8NUuj6ct9Y5K/VlWfSvK/k/xgd3+8qn61liZH+7kkP5zkZ6fbsBez9Hz5Wv5aktdV1Y1JPpXkZd39C1X1lUneNY0WP5zkm5N8dM6fJ0leneT10+/1kSQvP4PXntZKo+vd/Z6q+skkvzY1/Xh3vzdJqupNSd6XpdqPzrzsm5K8djrPn5OlZ7zfN6JGAJhH/eGdXAAAAMAT4XZxAAAAGETIBgAAgEGEbAAAABhEyAYAAIBBhGwAAAAYRMgGAACAQYRsAAAAGETIBgAAgEH+P2JEWn4ESy/qAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1224x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(17, 6))\n",
    "\n",
    "plt.hist(get_radial_distances(1000), alpha=0.3, density=True, rwidth=0.7)\n",
    "plt.hist(get_radial_distances(100), alpha=0.3, density=True, rwidth=0.7)\n",
    "plt.hist(get_radial_distances(10), alpha=0.3, density=True, rwidth=0.7)\n",
    "plt.hist(get_radial_distances(1), alpha=0.3, density=True, rwidth=0.7)\n",
    "\n",
    "plt.xlabel('Radial Distance from Mode')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bb38571",
   "metadata": {},
   "source": [
    "## Overthinking: Hamilton Monte Carlo in the raw"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e83b48f2",
   "metadata": {},
   "source": [
    "### R Code 9.5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e5cb3e4",
   "metadata": {},
   "source": [
    "### R Code 9.6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d6b2d43",
   "metadata": {},
   "source": [
    "### R Code 9.7"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25057b8a",
   "metadata": {},
   "source": [
    "### R Code 9.8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a7e9b7d",
   "metadata": {},
   "source": [
    "### R Code 9.9"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39cdc30c",
   "metadata": {},
   "source": [
    "### R Code 9.10"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbf37c09",
   "metadata": {},
   "source": [
    "## 9.4 Easy HMC: **ulam**"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}