Lecture15_03 - Sampling from a discrete distribution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Sampling from a discrete distribution\n",
    "\n",
    "##### Keywords: discrete sampling, MCMC, metropolis, poisson distribution, proposal matrix, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "hide": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "from scipy.stats import norm\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set_style(\"whitegrid\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In simulated annealing, we carried out combinatorical oprimization by sampling from a state space where each state was a vector of baseball simulation features.\n",
    "\n",
    "Since Metropolis MCMC is the same algorithm, it should be clear that its possible to simulate discrete possibilities in MCMC as long as you choose proposals which satisfy detailed balance. \n",
    "\n",
    "As an example, consider simulating a poisson distribution. Since its discrete, the proposal wont be a continuous $q(x,y)$ (the proposal probability to go from y to x), but rather a matrix indexed by a variable that corresponds to (indexes) the various states that can be obtained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def metropolis(p, qdraw, nsamp, xinit):\n",
    "    samples=np.empty(nsamp)\n",
    "    x_prev = xinit\n",
    "    accepted=0\n",
    "    for i in range(nsamp):\n",
    "        x_star = qdraw(x_prev)\n",
    "        p_star = p(x_star)\n",
    "        p_prev = p(x_prev)\n",
    "        pdfratio = p_star/p_prev\n",
    "        if np.random.uniform() < min(1, pdfratio):\n",
    "            samples[i] = x_star\n",
    "            x_prev = x_star\n",
    "            accepted+=1\n",
    "        else:#we always get a sample\n",
    "            samples[i]= x_prev\n",
    "            \n",
    "    return samples, accepted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Sampling a Poisson\n",
    "\n",
    "The poisson pmf is:\n",
    "\n",
    "$$p(k) = e^{-\\mu}\\frac{\\mu^k}{k!}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD3CAYAAADxJYRbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAZ30lEQVR4nO3dfZRc9X3f8beeWG2xHuJTHEXqAz2p8zmbBmRlZQq17GIszto8RFvSnrrUdoyPiHBSU1CiGOEH0tYttLLMU4x9rFjHNtShoZEUGQxCBUrQgko8iBORDF9nax/XXmEbUq+kot2xtFL/mDv4apjZvYPmce/n9dfM7/e7d793tPrs3N+9M785p06dwszM8mFupwswM7P2ceibmeWIQ9/MLEcc+mZmOeLQNzPLkfmdLmA6zz///Km+vr5OlzGtUqlEt9cIrrMVeqVW19lcvVDnsWPHXhkcHDynVl9Xh35fXx8DAwOdLmNaxWKx62sE19kKvVKr62yuXqizUCh8r16fp3fMzHLEoW9mliMOfTOzHHHom5nliEPfzCxHuvrunbzbdWCMLXuCQ+MTLF/az6YhMbxqxRns46U3tA8zmz1mDH1Jc4F7gJVACVgfEaNVY84BRoDzI2JS0k3Ae5PupcCyiFgm6UZgPfBy0rchIqI5hzK77DowxuYdB5k4PgXA2PgEm3ccBMgc2s3Yh5nNLlmmd4aBhRFxEXATsDXdKWkIeBRYVmmLiNsi4uKIuBj4AfChpGsQ+FClz4Ff35Y98VpYV0wcn2LLnuwvWTP2YWazS5bQXwM8AhAR+4HVVf0ngbXA/63eUNJVwE8i4tGkaRDYLGmfpM1vuOocODQ+0VB7q/ZhZrNLljn9xcDh1PMpSfMj4gRAROwFkFRr283Av0o9vx/4PHAE2Cnpioh4sN4PLpVKFIvFDCV2zuTkZEtqPOfs+fz41RM127P+vGbso91a9Xq2Qq/U6jqbq1fqrCdL6B8BFqWez60E/nQk/TIwXpn/lzQHuCMiDifPHwJWAXVDP89fw3DzFYtPm48H6F8wj5uv+BUGBrLNxzdjH+3WCx9xr+iVWl1nc/VCnYVCoW5fltAfAa4E/ljShcDBjD93LfBw6vli4AVJA8CrwCXA9oz7yp3KhdYzuXunGfsws9klS+jvBC6V9DQwB7hG0kZgNCJ2T7OdgL2VJxFxWNLNwBOU7wJ6LCK++cZLn/2GV60444Cu7KMX3p2YWevNGPoRcRK4rqr5xRrjzq16/ts1xtwL3NtYiWZm1iz+RK6ZWY449M3McsShb2aWIw59M7McceibmeWIQ9/MLEcc+mZmOeLQNzPLEYe+mVmOOPTNzHLEoW9mliMOfTOzHHHom5nliEPfzCxHHPpmZjni0DczyxGHvplZjjj0zcxyxKFvZpYjDn0zsxyZcWF0SXOBe4CVQAlYHxGjVWPOAUaA8yNiUtIc4AfAXydDnomIzZKuBD4NnAC2R8S25h2KmZnNZMbQB4aBhRFxkaQLga3AukqnpCHgNmBZaptfBJ6LiCtT4xYAtwNvB14FRiTtjogfnflhmJlZFllCfw3wCEBE7Je0uqr/JLAWKKTaBoEVkp4AJoAbgT5gNCJ+AiBpH/Au4IEzOgIzM8ssS+gvBg6nnk9Jmh8RJwAiYi+ApPQ2LwG3RsQDktYA91EO/vR+jgJLpvvBpVKJYrGYocTOmZyc7PoawXW2Qq/U6jqbq1fqrCdL6B8BFqWez60E/jS+RXnenojYJ2k55ZBP72cRMD7dTvr6+hgYGMhQYucUi8WurxHOrM5dB8bYsic4ND7B8qX9bBoSw6tWNLnCsl55PaF3anWdzdULdRYKhbp9We7eGQEuA0jm9A9m2OYW4IZkm5XA94G/At4q6c2SzqI8tfNMhn1ZB+06MMbmHQcZG5/gFDA2PsHmHQfZdWCs06WZ2RuQJfR3ApOSnqZ8IfZGSRsl/do029wG/FNJTwKfAz4cEceBjcAeymG/PSKcHF1uy55g4vjUaW0Tx6fYsic6VJGZnYkZp3ci4iRwXVXzizXGnZt6/BPg8hpjvgF8o+EqrWMOjU801G5m3c0fzrJpLV/a31C7mXU3h75Na9OQ6F8w77S2/gXz2DSkOluYWTfLcveO5VjlLp123b1jZq3l0G+Rdt7m2GrDq1b0bO1mdjqHfgtUbnOs3PVSuc0RcHiaWUd5Tr8FfJujmXUrh34L+DZHM+tWDv0W8G2OZtatHPot4Nsczaxb+UJuC/g2RzPrVg79FvFtjmbWjTy9Y2aWIw59M7McceibmeWIQ9/MLEcc+mZmOeLQNzPLEYe+mVmOOPTNzHLEoW9mliMzfiJX0lzgHmAlUALWR8Ro1ZhzgBHg/IiYlLQEuA9YDJwFbIyIZyT9M+CzwPeTTW+JiCebdjRmZjatLF/DMAwsjIiLJF0IbAXWVTolDQG3ActS22wEHouIOyQJ+CPgV4FB4Pci4k+adQBmZpZdltBfAzwCEBH7Ja2u6j8JrAUKqbbbKZ8VVH7GZPJ4EFgl6QbgWeDjEXHiDdZuZmYNyhL6i4HDqedTkuZXwjoi9gKU39CXRcR40raM8jTPDUnXXmAX8F3gi8B1wB/U+8GlUolisZj1WDpicnKy62sE19kKvVKr62yuXqmzniyhfwRYlHo+N8u7c0nnAfcDv5uat9+e+oPwp8CvT7ePvr4+BgYGMpTYOcVisetrBNfZCr1Sq+tsrl6os1Ao1O3LcvfOCHAZQDKnf3CmDST9MvAAcHVEPJy0zQH+QtLfSYa9h9OnhMzMrMWyvNPfCVwq6WlgDnCNpI3AaETsrrPNrcBC4M5k2udwRKyTtB7YIWkC+Ctg2xkfgZmZZTZj6EfEScpz72kv1hh3burxuur+pP1R4NHGSjQzs2bxh7PMzHLEoW9mliMOfTOzHHHom5nliEPfzCxHHPpmZjni0DczyxGHvplZjjj0zcxyxKFvZpYjDn0zsxxx6JuZ5YhD38wsRxz6ZmY54tA3M8uRLIuomJ2xXQfG2LInODQ+wfKl/WwaEsOrVnS6LLPccehby+06MMbmHQeZOD4FwNj4BJt3lFfddPCbtZend6zltuyJ1wK/YuL4FFv2RIcqMssvh7613KHxiYbazax1HPrWcsuX9jfUbmat49C3lts0JPoXzDutrX/BPDYNqUMVmeXXjBdyJc0F7gFWAiVgfUSMVo05BxgBzo+ISUn9wH3AW4CjwG9ExMuSrgQ+DZwAtkfEtqYejXWlysVa371j1nlZ7t4ZBhZGxEWSLgS2AusqnZKGgNuAZaltPgocjIjfl/R+4JOSfhe4HXg78CowIml3RPyoScdiXWx41QqHvFkXyBL6a4BHACJiv6TVVf0ngbVAoWqb/5I8fhj4FDAAjEbETwAk7QPeBTxQ7weXSiWKxWKGEjtncnKy62sE19kKvVKr62yuXqmzniyhvxg4nHo+JWl+RJwAiIi9AJLqbXMUWFJjP5X2uvr6+hgYGMhQYucUi8WurxFcZyv0Sq2us7l6oc5CoVC3L8uF3CPAovQ2lcDPuM0iYLzGfirtZmbWJllCfwS4DCCZ0z/YyDbA+4CngCLwVklvlnQW5amdZxqu2MzM3rAs0zs7gUslPQ3MAa6RtJHy/PzuOtt8AfhqMm//U+DqiDiebLeH8h+b7RExduaHYGZmWc0Y+hFxEriuqvnFGuPOTT0+BvyLGmO+AXyj4SrNzKwp/OEsM7McceibmeWIQ9/MLEcc+mZmOeLQNzPLEYe+mVmOOPTNzHLEoW9mliMOfTOzHMnyNQy5s+vAmBf8MLNZyaFfZdeBMTbvOMjE8SkAxsYn2Lyj/B1zDn4z63We3qmyZU+8FvgVE8en2LInOlSRmVnzOPSrHBqfaKjdzKyXOPSrLF/a31C7mVkvcehX2TQk+hfMO62tf8E8Ng2pzhZmZr3DF3KrVC7W+u4dM5uNHPo1DK9a4ZA3s1nJ0ztmZjni0DczyxGHvplZjsw4py9pLnAPsBIoAesjYjTVfy2wATgBfCYiHpR0B/C2ZMgyYDwiLpR0J7AGOJr0rYuIw007GjMzm1aWC7nDwMKIuEjShcBWYB2ApGXA9cBqYCGwT9LeiLgh6V8A7AOuTfY1CAxFxCvNPQwzM8siS+ivAR4BiIj9klan+i4ARiKiBJQkjQLnA3+e9H8MeDQiDiZnDG8FviTp54EvR8T26X5wqVSiWCw2dkRtNjk52fU1gutshV6p1XU2V6/UWU+W0F8MpKdgpiTNj4gTNfqOAksAJJ1FedrngqTvbOBu4HPAPOAJSd+KiL+o94P7+voYGBjIeiwdUSwWu75GcJ2t0Cu1us7m6oU6C4VC3b4sF3KPAIvS2ySBX6tvETCePF4L/Flqzv4YcGdEHIuIo8DjlK8TmJlZm2QJ/RHgMoBkTv9gqu9Z4J2SFkpaAgwALyR9a4GHU2N/CRiRNC+Z618DPHeG9ZuZWQOyhP5OYFLS08DtwI2SNkr6tYj4IXAX8BTld+6fiIjJZDsB36nsJCKKwL3AfuBJ4GsR8ZfNOxQzM5vJjHP6EXESuK6q+cVU/zZgW43tLq/RtgXY0niZZmbWDP5wlplZjjj0zcxyxN+yaT3DC9abnTmHvvUEL1hv1hye3rGe4AXrzZrDoW89wQvWmzWHQ996ghesN2sOh771BC9Yb9YcvpBrPcEL1ps1h0PfeoYXrDc7c57eMTPLEYe+mVmOOPTNzHLEoW9mliMOfTOzHHHom5nliEPfzCxHHPpmZjni0DczyxGHvplZjsz4NQyS5gL3ACuBErA+IkZT/dcCG4ATwGci4kFJbwa+DbyQDNsZEXfWGtvUozEzs2ll+e6dYWBhRFwk6UJgK7AOQNIy4HpgNbAQ2CdpL/CrwB9FxMcqO6k3NiJKzTwgMzOrL0vorwEeAYiI/ZJWp/ouAEaS4C5JGgXOBwaBQUlPAj+mHPZvrzP2z5t2NGZmNq0sob8YOJx6PiVpfkScqNF3FFgCvAgUIuJ/SPrXwN3Arjpj6yqVShSLxQwlds7k5GTX1wiusxV6pVbX2Vy9Umc9WUL/CLAo9XxuEvi1+hYB48D/Ao4lbTuBfw98rc7Yuvr6+hgYGMhQYucUi8WurxFcZyv0Sq2us7l6oc5CoVC3L8vdOyPAZQDJnP7BVN+zwDslLZS0BBigfPH2D4FfT8a8ByhMM9bMzNokyzv9ncClkp4G5gDXSNoIjEbEbkl3AU9R/gPyiYiYlHQTsF3SbwGvUr7j54e1xrbioMzMrLYZQz8iTgLXVTW/mOrfBmyr2ua7wLtr7Ot1Y83MrH384Swzsxxx6JuZ5YhD38wsRxz6ZmY54tA3M8sRh76ZWY449M3MciTLh7PMZo1dB8bYsic4ND7B8qX9bBoSw6tWdLoss7Zx6Ftu7DowxuYdB5k4PgXA2PgEm3eUv1XEwW954ekdy40te+K1wK+YOD7Flj3RoYrM2m/WvdP36bvVc2h8oqF2s9loVr3Tr5y+j41PcIqfnb7vOjDW6dKsCyxf2t9Qu9lsNKtC36fvNp1NQ6J/wbzT2voXzGPTkDpUkVn7zarpHZ++23Qq03ye/rM8m1Whv3xpP2M1At6n71YxvGqFQ95ybVZN7/j03cxserPqnb5P383MpjerQh98+m5mNp1ZNb1jZmbTc+ibmeXIjNM7kuYC9wArgRKwPiJGU/3XAhuAE8BnIuJBSX8P2J7sfw7wmxERkm4E1gMvJ5tviAjfRG9m1iZZ5vSHgYURcZGkC4GtwDoAScuA64HVwEJgn6S9wH8A/iAidkkaAm4FrgIGgQ9FRKH5h2JmZjPJMr2zBngEICL2Uw74iguAkYgoRcRhYBQ4H/gd4KFkzHxgMnk8CGyWtE/S5ibUb2ZmDcjyTn8xcDj1fErS/Ig4UaPvKLAkIl4BkCTgs5TPFgDuBz4PHAF2SroiIh6s94NLpRLFYjHzwXTC5ORk19cIrrMVeqVW19lcvVJnPVlC/wiwKPV8bhL4tfoWAeMAkt5N+VrAB5P5/DnAHckZAZIeAlYBdUO/r6+PgYGBrMfSEcVisetrBNfZCr1Sq+tsrl6os1CoP4OeJfRHgCuBP07m9A+m+p4F/qOkhUAfMAC8kAT+ncB7I+J7ydjFSd8A8CpwCeWLvWZm1iZZQn8ncKmkpynfiXONpI3AaETslnQX8BTl6wOfiIhJSXcAZwFfLc/wEBGxQdLNwBOU7wJ6LCK+2YJjMjOzOmYM/Yg4CVxX1fxiqn8bsK1qm5V19nUvcG/jZZp1j9MX6nnJX/VhPWXWfQ2DWSt5nV3rdf5ErlkDvFCP9TqHvlkDvFCP9TqHvlkDvM6u9TqHvlkDvFCP9TpfyDVrgBfqsV7n0DdrUGWhnl74ZKZZNU/vmJnliEPfzCxHHPpmZjni0DczyxFfyDXrgNO/v8d3AFn7OPTN2szf32Od5Okdszbz9/dYJzn0zdrM399jneTQN2szf3+PdZJD36zNmvH9PbsOjPGO2x7nH9z0EO+47XF2HRhrdpk2S/lCrlmbnen39/hCsJ0Jh75ZB1S+v+eNmO5CcCP79G2j+eTQN+sxzbgQ3KyzBf/h6D0zhr6kucA9wEqgBKyPiNFU/7XABuAE8JmIeFDS3wa+DvQDh4BrIuJYrbHNPiCz2W750n7GagR8IxeCm3G20Jo/HG9soflm/PGZTfuYTpYLucPAwoi4CLgJ2FrpkLQMuB54BzAE3CqpD/g08PWIeCdwANgwzVgza0AzLgQ342yhGZ83qPzhGBuf4BQ/+8PRyIVp76MxWUJ/DfAIQETsB1an+i4ARiKiFBGHgVHg/PQ2wMPA2mnGmlkDhlet4NarzmPF0n7mACuW9nPrVec19G6wGbeNdssfDu+jMVnm9BcDh1PPpyTNj4gTNfqOAkuq2mu1pdvrKpVKFIvFDCV2zuTkZNfXCK6zFTpZqxbCH677hVTLEYrFIzXH1qrz6vPexF1PT1KaOvVaW9+8OVx93psyH9M5Z8/nx6+eqNmedR/T/eHo1n3U+3fvlmOZSZbQPwIsSj2fmwR+rb5FwHiqfaJGW/XYuvr6+rp+ZaJeWT3JdTZfr9Raq86BAVix/Mzmjm++YvFpc/pQnma6+YpfYWAg236WL32p7vWJrK9tu/dR79+9W44FoFAo1O3LMr0zAlwGIOlC4GCq71ngnZIWSloCDAAvpLcB3gc8Nc1YM+uA4VUrGLnpEr572+WM3HRJwxcLmzHN1IzrE95HY7K8098JXCrpaWAOcI2kjcBoROyWdBflUJ8LfCIiJiV9BvhqcrfOK8DVEfFqrbFNOxIza7sz+bxBZXs4s4XmvY/GzDl16tTMozqkWCye6vbT514+xe9GvVIn9E6trrO5eqHOQqFQGBwcXF2rz9+9Y2aWIw59M7McceibmeWIQ9/MLEcc+mZmOdLVd+8UCoWXge91ug4zsx7z9wcHB8+p1dHVoW9mZs3l6R0zsxxx6JuZ5YhD38wsRxz6ZmY54tA3M8sRh76ZWY5k+Wrl3JO0ANgOnAv0UV7UfXeq/0ZgPfBy0rQhIpq3vlmDJD1HedEagO9GxDWpvq5YnF7Sh4EPJ08XAm8DlkXEeNJ/J+VlN48mY9Yly2y2s8Z/DPzniLhY0j8EvgKcorwOxG9HxMnU2H7gPuAtSc2/EREvv36vLa/zbcDdwBRQAj4UET+qGl/396PNta4CHgT+Oun+QkT8t9TYbnlN7weWJV3nAvsj4v2psXOAH/Cz43gmIja3o843wqGfzQeAv4mID0p6M/A8sDvVP0j5P1f95WraRNJCYE5EXFyjr7I4/WrKQbtP0t6IKLW3SoiIr1AOUSR9HtheCfzEIDAUEa+0u7akpt8DPgi8mjR9DvhkRPxPSV8E1lFea6Lio8DBiPh9Se8HPgn82w7UeSfwsYh4XtIG4OPAxtT4ur8fHah1EPhcRGyts0lXvKaVgJf0c8ATwI1Vm/wi8FxEXNnq2prB0zvZPAB8Knk8h/K75LRBYLOkfZI6/Rd+JfC3JD0q6fFktbOKrlucXtJq4B9FxJdSbXOBtwJfkjQi6SMdKO1/A1elng8CTyaPHwbWVo1fAzwyTX+rVNf5/oh4Pnk8H6heqGi6349Wq/WaXi7pzyR9WdKiqvHd8ppW/Dvg7oh4qap9EFgh6QlJ35TUvGWuWsChn0FE/L+IOJr8Uv53yu840u4HrgMuAdZIuqLdNaYcAz4LDCU1/VdJlTO6hhenb4ObKf9nSjub8hTFB4D3Ar8lqa1/nCLiT4DjqaY5EVH5+Hqt1y392rbtda2usxJIkv4J8G+A26s2me73o621Ul5CdVNEvAv4DnBL1SZd8ZoCSHoL8B6Ss9MqLwG3RsS7gf9EeUqqazn0M5L0dymf2t0bEV9Ptc8B7oiIVyLip8BDwKoOlQnwbeC+iDgVEd8G/gb4haSv4cXpW0nSUkAR8URV1zHgzog4FhFHgccpv0PtpJOpx7Vet/Rr2+nX9V8CXwQurzEHPt3vR7vtTE2J7uT1/2+65jUF/jnw9YiYqtH3LeBPASJiH7A8yYWu5NDPQNLPA48CH4+I7VXdi4EXJL0p+Ye+BOjk3P5HgK0AkpYn9VVOR7ttcfp3AY/VaP8lYETSvOQi+hrgubZW9noHJF2cPH4f5bWe00aAy6bpbwtJH6D8Dv/iiPhOjSHT/X602x5JFySP38Pr/990xWuaWEt5iqmWW4AbACStBL6fOivsOr6Qm83NwM8Bn5JUmdvfBpwdEV+SdDPls4AS8FhEfLNDdQJ8GfiKpH2U7zT5CHC9pLoL2XewVlE+rS8/kTYClTrvBfZTPs3+WkT8ZYdqrPgdYJuks4Ai5Wk+JD0KXAF8Afhq8rr/FLi63QVKmgfcBfwfYEcytfxkRNwi6WuUpyVf9/sREdXXqNrlo8Ddko4DPwR+E7rrNU057XcVTqvzNuA+SZdTvt734bZX1wB/y6aZWY54esfMLEcc+mZmOeLQNzPLEYe+mVmOOPTNzHLEoW9mliMOfTOzHPn/wImYB2oThAcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.stats import poisson\n",
    "xxx= np.arange(1,20,1)\n",
    "plt.plot(xxx, poisson.pmf(xxx, mu=5), 'o'); "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To sample from this distribution, we must create a proposal matrix which allows us to go from any integer output to any other in a finite number of steps. This matrix must be symmetric, since we wish to use Metropolis.\n",
    "\n",
    "A simple such matrix, which is although a bit slow, would be one which has immediate off-diagonal elements (from Stats 580 at Iowa state..)\n",
    "\n",
    "![](images/propmatrix.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prop_pdf(ito, ifrom):\n",
    "    if ito == ifrom - 1:\n",
    "        return 0.5\n",
    "    elif ito == ifrom + 1:\n",
    "        return 0.5\n",
    "    elif ito == ifrom and ito == 0:#needed to make first row sum to 1\n",
    "        return 0.5\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prop_draw(ifrom):\n",
    "    u = np.random.uniform()\n",
    "    if ifrom !=0:\n",
    "        if u < 1/2:\n",
    "            ito = ifrom -1\n",
    "        else:\n",
    "            ito = ifrom + 1\n",
    "    else:\n",
    "        if u < 1/2:\n",
    "            ito=0\n",
    "        else:\n",
    "            ito=1\n",
    "    return ito"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "rv = poisson(5)\n",
    "samps, acc = metropolis(rv.pmf, prop_draw, 50000, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "41389"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "acc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD3CAYAAADxJYRbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAcQElEQVR4nO3df3Rc5X3n8ffIkiXD2nIxZL2uk/qQTb9Rl9ihcli8R8mSGA4pxJhNd73ZbJKGHFOTpvGC2yWGpMDStDhuiPnRkJyY+iTAOkQpdSJDILgOS7HASyPCkWnH31QhpHadpEBj2Ys8Y8nS/nHvhKvRjHTHujOj0f28/tHc53nu1feOrj5z55mZO5mxsTFERCQdmupdgIiI1I5CX0QkRRT6IiIpotAXEUkRhb6ISIo017uAyTz//PNjra2tVdt+Pp+nmttPgmpMhmpMhmpMTjXrHBoaeqWzs/OcUn0zOvRbW1vp6Oio2vaz2WxVt58E1ZgM1ZgM1ZicatbZ19f3k3J9mt4REUkRhb6ISIoo9EVEUkShLyKSIgp9EZEUUeiLiKSIQl9EJEUU+iIiKaLQFxFJEYX+LJcbPlW6o78btp0HtywMfvZ3n/Y23vqNVVNuQ0Rmhikvw2BmTcA9wAogD6x394GiMecAvcByd8+Z2WbgvWH3QmCxuy82s+uA9cDLYd8Gd/dkdkVKaWuZw7LNj4xru6JpH1ta7uWMzMmgYfAQQw99gs0P/oCe0a4J23hpy+UTN9zfDbs3wvAJMuE22L0x6Fu+LtmdEJHExDnTvxJoc/dVwGbg9minmV0KPA4sLrS5+xZ3v8jdLwIOAx8JuzqBjxT6FPj1cX1z9+uBHzojc5Lrmys4U997KwyfGN82fCJoF5EZK84F17qAxwDcfb+ZrSzqHwUuBvqKVzSz9wO/cPfHw6ZO4AYzWww84u63TfaL8/k82Ww2RomnJ5fLVXX7SZhujaUu6LQk80rJsUsyr5bdTnENbx08HJzhFxkbPMzBGXifpuFvXQuqMTn1qjNO6C8ABiPLp8ys2d1HANx9D4CZlVr3BuC/RZYfBL4IHAN2mdn73P3hcr9YV9msTo1Hxs5maYngPzK2qOT43PCpiTW0Lw2mdIpk2pdOGJsbPkVby5zTLzgBaf1bJ001JqfKV9ks2xcn9I8B8yPLTYXAn4yZ/QZwtDD/b2YZ4A53HwyXHwHOB8qGvlTH1pF14+f0gaGxuWwdKT0XX/p1gTUlt7H55TX0FI0t+ZqAiNRFnNDvBdYA3WZ2IXAg5rYvBh6NLC8AXjCzDuA14D3AjgpqlYT0jHbBcDC3vyTzKkfGFrF1ZF3JF3GruQ0Rqb04ob8LuMTMngYywFVmtgkYcPeeSdYzYE9hwd0HzexG4AmCdwHtdffvnH7pMh09o130nJxeQCexDRGprSlD391HgWuKmg+WGLesaPkTJcbcD9xfWYkiIpIUfThLRCRFFPoiIimi0BcRSRGFvohIiij0RURSRKEvIpIiCn0RkRRR6IuIpIhCX0QkRRT6IiIpotAXEUkRhb6ISIoo9EVEUkShLyKSIgp9EZEUUeiLiKSIQl9EJEUU+iIiKaLQFxFJEYW+iEiKKPRFRFKkeaoBZtYE3AOsAPLAencfKBpzDtALLHf3nJllgMPAP4RDnnH3G8xsDXATMALscPftye2KiIhMZcrQB64E2tx9lZldCNwOrC10mtmlwBZgcWSdNwPPufuayLgWYBvwDuA1oNfMetz959PfDRERiSNO6HcBjwG4+34zW1nUPwpcDPRF2jqBXzWzJ4ATwHVAKzDg7r8AMLN9wLuAb5b7xfl8nmw2G3NXKpfL5aq6/SRMt8aOjo4Eqzl99b6f0/C3rgXVmJx61Rkn9BcAg5HlU2bW7O4jAO6+B8DMouv8FLjN3b9pZl3AAwTBH93OcaB9sl/c2tpa1dDKZrMzJhTLaYQa46j3PjTC/agak9EINUJ16+zr6yvbF+eF3GPA/Og6hcCfxPeBbwO4+z5gCUHIR7czHzga4/fLbNXfDdvOg1sWBj/7u+tdkcisFyf0e4HLAMI5/QMx1rkZuDZcZwVwCPh74C1mdpaZzSWY2nnmdIqWWaC/G3ZvhMFDwFjwc/dGBb9IlcUJ/V1AzsyeJngh9joz22RmV0yyzhbgP5rZk8AXgI+6+zCwCfguQdjvcPd/ml750rD23grDJ8a3DZ8I2kWkaqac03f3UeCaouaDJcYti9z+BXB5iTG7gd0VVymzz+DhytpFJBH6cJbUR/vSytpFJBEKfamP1TdBy7zxbS3zgnYRqRqFvlRdbvjUxMbl62DNXdD+RiAT/FxzV9AedxsiUrE479OXmai/O3jRc/BwMCWy+qaygVlvbS1zWLb5kRI9ZwKfC27mgJ3AzlLj4KUtE14iEpHToNCf4d607NyJjYW3Oxbe/VJ4uyPM2OAXkZlBoT/DnTmvdcJZ8r65N7C0aeLbHQ//5Q107TxzXLPOkEUkSnP6DWhJ5pUy7a/WuBIRaTQK/QZ0ZOzsMu2LalyJiDQahX4D2jqyjqGxuePahsbmsnVE8/kiMjnN6TegntEuGIbrm7tZknmVI2OL2DqyLmgXEZmEQr9B9Yx20XNSIS8ildH0johIiij0RURSRKEvIpIiCn0RkRRR6IuIpIhCX0QkRRT6IiIpotAXEUkRhb6ISIoo9EVEUmTKyzCYWRNwD7ACyAPr3X2gaMw5QC+w3N1zZtYOPAAsAOYCm9z9GTP7T8DngUPhqje7+5OJ7Y2IiEwqzrV3rgTa3H2VmV0I3A6sLXSa2aXAFmBxZJ1NwF53v8PMDPg68JtAJ3C9uz+U1A6IiEh8cUK/C3gMwN33m9nKov5R4GKgL9K2jeBZQeF35MLbncD5ZnYt8CzwKXcfOc3aRUSkQnFCfwEwGFk+ZWbNhbB29z0AwQl9wN2Phm2LCaZ5rg279gDfAn4MfBm4Bvjzcr84n8+TzWbj7kvFcrlcVbefhI6OjnqXMGNM52/VCH9r1ZiMRqgR6ldnnNA/BsyPLDfFOTs3s7cBDwJ/GJm33xF5QPg28NuTbaO1tbWqoZfNZhWqDWQ6f6tG+FurxmQ0Qo1Q3Tr7+vrK9sV5904vcBlAOKd/YKoVzOw3gG8CH3T3R8O2DNBvZkvDYasZPyUkIiJVFudMfxdwiZk9DWSAq8xsEzDg7j1l1rkNaAPuDKd9Bt19rZmtB/7KzE4Afw9sn/YeiIhIbFOGvruPEsy9Rx0sMW5Z5Pba4v6w/XHg8cpKFBGRpOjDWSIiKaLQFxFJEYW+iEiKKPRFRFJEoS8ikiIKfRGRFFHoi4ikiEJfRCRFFPoiIimi0BcRSRGFvohIiij0RURSRKEvIpIiCn0RkRRR6IuIpIhCXxpbfzdsOw9uWRj87O+ud0UiM1qcb84SmZn6u2H3Rhg+ESwPHgqWAZavq19dIjOYzvSlce299fXALxg+EbSLSEkKfWlcg4craxcRhb40sPallbWLiEJfGtjqm6Bl3vi2lnlBu4iUNOULuWbWBNwDrADywHp3Hygacw7QCyx395yZzQMeAN4AHAd+x91fNrM1wE3ACLDD3bcnujcya+WGT9HWMmd8Y+HF2r23BlM67UuDwC/xIu6blp1bgypFZr447965Emhz91VmdiFwO7C20GlmlwJbgMWRdT4OHHD3W8zsA8BnzOwPgW3AO4DXgF4z63H3nye0LzKLtbXMYdnmR0r0nAl8LriZA3YCOyeOe2nL5VWsTqRxxJne6QIeA3D3/cDKov5R4GLgX0qtAzwa9ncAA+7+C3c/CewD3nX6pYuISKXinOkvAAYjy6fMrNndRwDcfQ+AmZVb5zjQXmI7hfay8vk82Ww2RomnJ5fLVXX7Sejo6Kh3CbPGTP9bN8LxqBqTU68644T+MWB+ZLmpEPgx15kPHC2xnUJ7Wa2trVUNvWw2q1BNkZn+t26E41E1Jqeadfb19ZXtizO90wtcBhDO6R+oZB3gt4CngCzwFjM7y8zmEkztPBNjWyIikpA4Z/q7gEvM7GkgA1xlZpsI5ud7yqzzJeBrZrYPOAl80N2Hw/W+S/Bgs8Pd/2n6uyAiInFNGfruPgpcU9R8sMS4ZZHbQ8B/KTFmN7C74ipFRCQR+nCWiEiKKPRFRFJEoS8ikiIKfRGRFFHoi4ikiEJfRCRFFPoiIimi0BcRSRGFvohIiij0RURSRKFfL/3dsO08uGVh8LO/u94ViUgKxLngmiStvxt2b4ThE8Hy4KFgGUp+1Z+ISFJ0pl9FueFTpTv23vp64BcMnwjaRUSqSGf6VVTue11fbD1EU2bi+NGjhzm3aLy+21VEkqQz/To4MnZ2mfZFNa5ERNJGoV8HW0fWMTQ2d1zb0Nhcto5oPl9EqkvTO3XQM9oFw3B9czdLMq9yZGwRW0fWBe0iIlWk0K+TntEuek4q5EWktjS9IyKSIgp9EZEUUeiLiKTIlHP6ZtYE3AOsAPLAencfiPRfDWwARoDPuvvDZnYH8PZwyGLgqLtfaGZ3Al3A8bBvrbsPJrY3IiIyqTgv5F4JtLn7KjO7ELgdWAtgZouBjcBKoA3YZ2Z73P3asL8F2AdcHW6rE7jU3V9JdjdERCSOONM7XcBjAO6+nyDgCy4Aet09H56xDwDLI/2fBB539wPhM4a3AF8xs14z+1gieyAiIrHFOdNfAESnYE6ZWbO7j5ToOw60A5jZXIJpnwvCvjOBu4EvAHOAJ8zs++7eX+4X5/N5stls3H2pWC6Xq+r2Ozo6qrZtqVw1/9ZJqPbxmATVmJx61Rkn9I8B8yPLTWHgl+qbDxwNb18M/E1kzn4IuNPdhwDM7HsErxOUDf3W1taqBmc2m1Uwp8hM/1s3wvGoGpNTzTr7+vrK9sWZ3ukFLgMI5/QPRPqeBd5pZm1m1g50AC+EfRcDj0bG/jrQa2Zzwrn+LuC5uDshIiLTFyf0dwE5M3sa2AZcZ2abzOwKd/8ZcBfwFPA94NPungvXM+DFwkbcPQvcD+wHngTuc/e/S25XRERkKlNO77j7KHBNUfPBSP92YHuJ9SZcE9jd/wz4s8rLFBGRJOjDWSIiKaLQFxFJEYW+iEiKKPQl3fq7Ydt5cMvC4Gd/d70rEqkqXU9f0qu/G3ZvfP1L6gcPBcsAy/UtZjI76Uxf0mvvra8HfsHwiaBdZJZS6Et6DR6urF1kFlDoS3q1L62sXWQWUOhLeq2+CVrmjW9rmRe0i8xSCn1JhdzwqYmNy9fBmrug/Y1AJvi55q6SL+KWXF+kAendO5IKbS1zWLb5kRI9ZwKfC27mgJ3AzonjXtoy4aoiIg1JZ/oiIimi0BcRSRGFvohIiij0RURSRKEvIpIiCn0RkRRR6IuIpIhCX0QkRRT6IiIpotAXEUmRKS/DYGZNwD3ACiAPrHf3gUj/1cAGYAT4rLs/bGZnAT8EXgiH7XL3O0uNTXRvRERkUnGuvXMl0Obuq8zsQuB2YC2AmS0GNgIrgTZgn5ntAX4T+Lq7f7KwkXJj3T2f5A6JiEh5cUK/C3gMwN33m9nKSN8FQG8Y3HkzGwCWA51Ap5k9CfwzQdi/o8zYv01sb0REZFJxQn8BMBhZPmVmze4+UqLvONAOHAT63P2vzey/A3cD3yoztqx8Pk82m41R4unJ5XJV3X5HR0fVti21V81jBap/PCZBNSanXnXGCf1jwPzIclMY+KX65gNHgf8LDIVtu4BbgfvKjC2rtbW1qsGZzWYVzBJbtY+VRjgeVWNyqllnX19f2b44797pBS4DCOf0D0T6ngXeaWZtZtYOdBC8eHsv8NvhmNVA3yRjRUSkRuKc6e8CLjGzp4EMcJWZbQIG3L3HzO4CniJ4APm0u+fMbDOww8x+D3iN4B0/Pys1tho7JSIipU0Z+u4+ClxT1Hww0r8d2F60zo+Bd5fY1oSxIiJSO/pwlohIiij0RURSRKEvIpIiCn0RkRRR6IuIpIhCX0QkRRT6IiIpotAXEUkRhb7IdPV3w7bz4JaFwc/+7npXJFJWnMswiEg5/d2weyMMnwiWBw8FywDL19WvLpEydKYvMh17b3098AuGTwTtIjOQQn8SueFTpTv0dD51yh4Lg4djt5fdhkgNaXpnEm0tc1i2+ZFxbVc07WNLy72ckTkZNAweYuihT7D5wR/QM9o1buxLWy6vValSZaWOBYB9cxextOmVCe2HRxfRVTRex4PMBDrTr9D1zd2vB37ojMxJrm/W2X4abR1Zx9DY3HFtQ2Nz2Tqi+XyZmXSmX6ElmYlndUH7qzWuRGaCntEuGA5OBpZkXuXI2CK2jqyb8KxPZKZQ6FfoyNjZLC0R/EfGFtWhGpkJeka76DmpkJfGoOmdCunpvIg0Mp3pV0hP50WkkSn0T4OezotIo9L0johIiij0RURSZMrpHTNrAu4BVgB5YL27D0T6rwY2ACPAZ939YTN7E7Aj3H4G+F13dzO7DlgPvByuvsHdPckdEhGR8uLM6V8JtLn7KjO7ELgdWAtgZouBjcBKoA3YZ2Z7gD8G/tzdv2VmlwK3Ae8HOoGPuHtf8rsiIiJTiTO90wU8BuDu+wkCvuACoNfd8+4+CAwAy4E/AAqfQW8GcuHtTuAGM9tnZjckUL+IiFQgzpn+AmAwsnzKzJrdfaRE33Gg3d1fATAzAz5P8GwB4EHgi8AxYJeZvc/dHy73i/P5PNlsNvbOVCqXy026/Y6Ojqr9bkmnyY63qY7HmUA1JqdedcYJ/WPA/MhyUxj4pfrmA0cBzOzdBK8FfDicz88Ad4TPCDCzR4DzgbKh39raWtXgzWazCnapqcmOt0Y4HlVjcqpZZ19f+Rn0OKHfC6wBusM5/QORvmeBPzGzNqAV6ABeCAP/TuC97v6TcOyCsK8DeA14D8GLvSIiUiNxQn8XcImZPU3wTpyrzGwTMODuPWZ2F/AUwesDn3b3nJndAcwFvhbM8ODuvsHMbgSeIHgX0F53/04V9klERMqYMvTdfRS4pqj5YKR/O7C9aJ0VZbZ1P3B/5WWKiEgS9OEskXoLv4ntrd9YpW9ik6rTtXdE6inyxeoZ0BerS9XpTF+knvTF6lJjCn2RGin5xej6YnWpMU3viNRIqS9X1xerS63pTF+kjvRNbFJrOtMXqSN9E5vUmkJfpM70TWxSS5reERFJEYW+iEiKKPRFRFJEoS8yG4SXcuCWhbqUg0xKL+SKNLrIpRwAXcpBJpXOM31d4EpmE13KQSqQvjN9XeBKZpsKLuUgMmvP9Mtep0RnRdKgyh7T7Utjt+v6PTJrz/RLXecE4MXWQzRlJo4fPXqYc3WtE5nByh3TVzStYUvLvZyROfnLtqGxuWx+eQ09OqalyKw90y/nyNjZZdoX1bgSkWT0jHaxeXg9h0fPZnQsw+HRs9k8vL6ySzno3T+pMWvP9MvZOrKu5FmRLnAljWxal3LQu39SJXVn+omcFYk0qJJz+hW8zvWmZeeW37ieLTSE1J3pgy5wJelV6nWBRF7nSuLZQn938EAzeDh4EXr1TXqmUQVThr6ZNQH3ACuAPLDe3Qci/VcDG4AR4LPu/rCZnQ3sBOYBR4Cr3H2o1Nikd0hEKnNk7GyWZiZ+kUup17lyw6doa5kzcSOTPVsoCu6S20jwQeOt03nQmO4DTxIPXFV+8Itzpn8l0Obuq8zsQuB2YC2AmS0GNgIrgTZgn5ntAW4Cdrr7V81sM7DBzL5eaqy75xPbGxGpWCWvcyXxrriDf/zeiQMTfNA47c/fTPeBJ6kHriq/vhIn9LuAxwDcfb+ZrYz0XQD0hsGdN7MBYHm4zp+GYx4Nb/+ozNi/TWRPROS0JPFFLpU8W5juFNN0HzSgzANHhduYbg1V28YUMmNjY5MOMLN7gYfc/dFw+R+Bc919xMw+BLzN3T8V9t0H3Ad8OWw/YWbnFrWNG+vuf13ud/f19b0M/GS6Oyki1dWcP3rW3BP//GswFnlzSGb05Lw3/GSkdeG/TLX+vMEfvS0zNjK3uH0s03zyRPubD0y1/hlHvbNc39BC65tq/SS2MRNqiPi1zs7Oc0p1xDnTPwbMjyw3uftImb75wNFI+4kSbcVjyypXtIjMNmWzrkbrz5YaphbnLZu9wGUA4Zx+9FH3WeCdZtZmZu1AB/BCdB3gt4CnJhkrIiI1Emd6p/DuneVABriKINAH3L0nfEfO7xI8gPypuz9kZv8a+BrB2fwrwAfd/bVSY6u0XyIiUsKUoS8iIrNH6j6RKyKSZgp9EZEUUeiLiKTIrL/2zulcRqIONbYAO4BlQGtYR0+k/zpgPfBy2LTB3b0OdT5H8NZbgB+7+1WRvrrfj2EdHwU+Gi62AW8HFrv70bD/ToIPDx4Px6x198Ea1vfvgc+5+0Vm9m+BrwJjBO9k+4S7j0bGzgMeAN4Q1vs77v7yxK1Wtca3A3cDpwj+fz7i7j8vGl/2uKhRjecDDwP/EHZ/yd2/ERk7E+7HB4HFYdcyYL+7fyAyNgMcjuzDM+5+QzXqmvWhz2lcRqIOl4b4EPCqu3/YzM4Cngd6Iv2dBP9slXw4I1Fm1gZk3P2iEn0z5X7E3b9KEKSY2ReBHYXAD3UCl7r7xI+PVpmZXQ98GHgtbPoC8Bl3/z9m9mWC43JXZJWPAwfc/RYz+wDwGeB/1LjGO4FPuvvzZrYB+BSwKTK+7HFRwxo7gS+4++1lVqn7/VgIeDP7FeAJ4LqiVd4MPOfua6pZF6RjemfcZSQIgqngl5eRCM/2CpeGqLVvAn8U3s4QnC1HdQI3mNk+M6vKo38MK4AzzOxxM/te+ABaMFPux18KLxfy79z9K5G2JuAtwFfMrNfMPlbjsn4EvD+y3Ak8Gd5+FLi4aPwvj90y/dVQXOMH3P358HYzkCsaP9lxUasaO4HLzexvzOwvzGx+0fiZcD8W/C/gbnf/aVF7J/CrZvaEmX3HzKxahaUh9BcA0afvp8ysuUzfcaC9VoUVuPv/c/fj4cH6lwRnIlEPAtcA7wG6zOx9ta4RGAI+D1wa1vK/Z9r9WORGgn+wqDMJpio+BLwX+D0zq9mDU/i5lOFIU8bdC++ZLnWfRe/XmtynxTUWwsnM/gPw+8C2olUmOy5qUiPBBz//p7u/C3gRuLlolbrfjwBm9gZgNeEz0SI/BW5z93cTXKvsgWrVlobQP53LSNScmb2R4Gnf/e6+M9KeAe5w91fc/STwCHB+HUr8IfCAu4+5+w+BV4F/E/bNmPsRwMwWAubuTxR1DQF3uvuQux8Hvkdwplovo5Hbpe6z6P1az2PzvxJcO+vyEnPhkx0XtbIrMvW5i4n/HzPifgT+M8HVh0t9O/33gW8DuPs+YEn4v5+4NIT+6VxGoqbCTzA/DnzK3XcUdS8AXjCzfxUeBO8B6jG3/zGC10MwsyVhXYWnqDPifox4F7C3RPuvA71mNid88bwLeK6mlY33AzO7KLxduFxJVKnLmdRUeFHF3wcucvcXSwyZ7Liole+a2QXh7dVM/P+o+/0YuphgeqmUm4FrAcxsBXAo8iwwUWl4IXcXcImZPU14GQkz28Trl5G4i+AgaAI+7e7Fc5a1cCPwK8AfmVlhbn87cKa7f8XMbiR4FpAH9rr7d+pQ418AXzWzfQTvNvkYsNHMZtL9WGAET/ODhfF/7/uB/QRPve9z97+rU40AfwBsN7O5QJZgag8zexx4H/Al4GvhfX4S+GAtizOzOcBdwD8CfxVOMz/p7jeHV8n9DCWOi8gz6Vr5OHC3mQ0DPyO41MuMuR8jxh2XMK7GLcADZnY5wWt6H61WEboMg4hIiqRhekdEREIKfRGRFFHoi4ikiEJfRCRFFPoiIimi0BcRSRGFvohIivx/ErD66ufKKhIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xxx = np.arange(0,samps.max())\n",
    "plt.hist(samps, bins=xxx, density=True, align='left');\n",
    "plt.plot(xxx, rv.pmf(xxx),'o');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Raw Cell Format",
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}