{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using `numpy.random`\n",
    "\n",
    "- [Reference](https://docs.scipy.org/doc/numpy/reference/routines.random.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random number generators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Continuous distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Uniform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.1824764 ,  0.54965542,  0.47202485,  0.95005208],\n",
       "       [ 0.94097831,  0.02402875,  0.63758556,  0.37598786],\n",
       "       [ 0.15309177,  0.32923875,  0.36835565,  0.71566394]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.rand(3,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.79982921,  0.49167894,  0.31763543,  0.80053871],\n",
       "       [ 0.73086786,  0.82912843,  0.09356568,  0.25450986],\n",
       "       [ 0.75515315,  0.06455766,  0.10129754,  0.08602764]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.random_sample((3,4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.94293588,  0.51498122,  0.38971777,  0.59287419],\n",
       "       [ 0.18248637,  0.02499295,  0.842292  ,  0.57454183],\n",
       "       [ 0.81824035,  0.27704791,  0.62583653,  0.50282242]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.uniform(0, 1, (3,4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Normal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.60069026,  0.09858881, -1.4709185 , -0.70565061],\n",
       "       [-0.46729956, -1.6822813 , -0.29514808,  1.80174322],\n",
       "       [ 2.231532  ,  0.63953851, -0.78135108,  1.13137354]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.randn(3, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.0417215 ,  1.02702489,  0.88650073,  0.3359106 ],\n",
       "       [ 0.37837503,  1.37917244,  1.26210605, -1.01897335],\n",
       "       [ 0.65942071,  0.24987618,  0.05533211, -1.13927286]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.standard_normal((3,4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.1300124 , -1.23111144,  0.82896809,  0.17212567],\n",
       "       [-0.98791995, -0.23104989,  2.19710114, -1.26945484],\n",
       "       [-0.69139685,  0.50303272,  1.17180097, -0.54426858]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.normal(0, 1, (3,4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discrete distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Uniform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[6, 7, 7, 2],\n",
       "       [9, 6, 3, 5],\n",
       "       [2, 0, 4, 4]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.randint(0, 10, (3,4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Binomial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[3, 1, 3, 5],\n",
       "       [5, 2, 2, 3],\n",
       "       [3, 1, 1, 1]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.binomial(10, 0.3, (3,4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Poisson"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 0, 0, 0],\n",
       "       [0, 0, 0, 0],\n",
       "       [0, 2, 0, 1]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.poisson(0.2, (3,4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multivariate distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multinomial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have 3 fair 6-sided dice. You roll each of them 10 times. How many times do the faces 1,2,3,4,5,6 appear?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 3, 2, 1, 1, 2],\n",
       "       [0, 1, 2, 2, 3, 2],\n",
       "       [1, 2, 3, 0, 2, 2]])"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.multinomial(10, np.ones(6)/6, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multivariate normal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.array([[1, 2]]).T\n",
    "sigma = x @ x.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu = np.array([5,5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2],\n",
       "       [2, 4]])"
      ]
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 6.68175511,  8.3635102 ],\n",
       "       [ 4.11165655,  3.22331307],\n",
       "       [ 4.67467303,  4.3493461 ]])"
      ]
     },
     "execution_count": 115,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.multivariate_normal(mu, sigma, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Seed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.0856306 ,  0.99734545,  0.2829785 , -1.50629471],\n",
       "       [-0.57860025,  1.65143654, -2.42667924, -0.42891263],\n",
       "       [ 1.26593626, -0.8667404 , -0.67888615, -0.09470897]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.seed(123)\n",
    "np.random.randn(3,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.0856306 ,  0.99734545,  0.2829785 , -1.50629471],\n",
       "       [-0.57860025,  1.65143654, -2.42667924, -0.42891263],\n",
       "       [ 1.26593626, -0.8667404 , -0.67888615, -0.09470897]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.seed(123)\n",
    "np.random.randn(3,4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RNG class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rng1 = np.random.RandomState(123)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.0856306 ,  0.99734545,  0.2829785 , -1.50629471],\n",
       "       [-0.57860025,  1.65143654, -2.42667924, -0.42891263],\n",
       "       [ 1.26593626, -0.8667404 , -0.67888615, -0.09470897]])"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rng1.normal(0, 1, (3,4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.49138963, -0.638902  , -0.44398196, -0.43435128],\n",
       "       [ 2.20593008,  2.18678609,  1.0040539 ,  0.3861864 ],\n",
       "       [ 0.73736858,  1.49073203, -0.93583387,  1.17582904]])"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rng1.normal(size=(3,4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rng2 = np.random.RandomState(123)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.0856306 ,  0.99734545,  0.2829785 , -1.50629471],\n",
       "       [-0.57860025,  1.65143654, -2.42667924, -0.42891263],\n",
       "       [ 1.26593626, -0.8667404 , -0.67888615, -0.09470897]])"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rng2.normal(0, 1, (3,4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.49138963, -0.638902  , -0.44398196, -0.43435128],\n",
       "       [ 2.20593008,  2.18678609,  1.0040539 ,  0.3861864 ],\n",
       "       [ 0.73736858,  1.49073203, -0.93583387,  1.17582904]])"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rng2.normal(size=(3,4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 1, 3, 2, 0, 2, 2, 3, 0, 3])"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.choice(4, size=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([7, 5, 7, 3, 7, 5, 3, 7, 3, 5])"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.choice([1,3,5,7], 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['T', 'T', 'T', 'H', 'T', 'T', 'T', 'T', 'H', 'H'],\n",
       "      dtype='<U1')"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.choice(['H', 'T'], 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['A', 'G', 'T', 'C', 'A', 'G', 'T', 'T', 'T', 'T'],\n",
       "      dtype='<U1')"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.choice(list('ACTG'), 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'CTAGGGCATA'"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "''.join(np.random.choice(list('ACTG'), 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.arange(24).reshape((6,4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3],\n",
       "       [ 4,  5,  6,  7],\n",
       "       [ 8,  9, 10, 11],\n",
       "       [12, 13, 14, 15],\n",
       "       [16, 17, 18, 19],\n",
       "       [20, 21, 22, 23]])"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[20, 21, 22, 23],\n",
       "       [16, 17, 18, 19],\n",
       "       [ 4,  5,  6,  7],\n",
       "       [16, 17, 18, 19],\n",
       "       [16, 17, 18, 19],\n",
       "       [12, 13, 14, 15],\n",
       "       [16, 17, 18, 19],\n",
       "       [12, 13, 14, 15],\n",
       "       [ 0,  1,  2,  3],\n",
       "       [20, 21, 22, 23]])"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ridx = np.random.choice(x.shape[0], 10)\n",
    "x[ridx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  0,  3,  0,  0,  3,  0,  0,  1,  0],\n",
       "       [ 4,  4,  7,  4,  4,  7,  4,  4,  5,  4],\n",
       "       [ 8,  8, 11,  8,  8, 11,  8,  8,  9,  8],\n",
       "       [12, 12, 15, 12, 12, 15, 12, 12, 13, 12],\n",
       "       [16, 16, 19, 16, 16, 19, 16, 16, 17, 16],\n",
       "       [20, 20, 23, 20, 20, 23, 20, 20, 21, 20]])"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cidx = np.random.choice(x.shape[1], 10)\n",
    "x[:, cidx]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Permutations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In-place"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.arange(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.random.shuffle(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([9, 1, 3, 6, 0, 5, 7, 2, 4, 8])"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make a copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.arange(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 3, 9, 1, 5, 4, 2, 7, 8, 6])"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.permutation(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Can also give the size to get permutation indexes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 0, 6, 8, 9, 3, 4, 5, 7, 1])"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.permutation(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1**. Simulate 100 tosses of a fair coin. What is the longest run of consecutive heads and where in the sequence does it occur?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tosses = np.random.choice(['H', 'T'], 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['T', 'H', 'H', 'H', 'T', 'H', 'T', 'T', 'T', 'H', 'H', 'T', 'H',\n",
       "       'T', 'H', 'T', 'T', 'T', 'T', 'H', 'T', 'H', 'H', 'H', 'T', 'H',\n",
       "       'T', 'H', 'T', 'T', 'H', 'H', 'T', 'H', 'H', 'T', 'T', 'T', 'T',\n",
       "       'T', 'H', 'H', 'T', 'T', 'T', 'H', 'T', 'H', 'T', 'T', 'H', 'H',\n",
       "       'H', 'H', 'H', 'T', 'H', 'H', 'T', 'H', 'H', 'H', 'H', 'T', 'T',\n",
       "       'T', 'T', 'H', 'H', 'H', 'T', 'T', 'H', 'T', 'H', 'H', 'H', 'H',\n",
       "       'T', 'H', 'H', 'H', 'H', 'H', 'T', 'T', 'H', 'H', 'H', 'H', 'T',\n",
       "       'T', 'H', 'T', 'H', 'H', 'H', 'T', 'H', 'T'],\n",
       "      dtype='<U1')"
      ]
     },
     "execution_count": 147,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tosses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def find_runs(xs, target='H'):\n",
    "    count = 0\n",
    "    max_count = 0\n",
    "    idx = 0\n",
    "    for i, x in enumerate(xs):\n",
    "        if x == target:\n",
    "            count += 1\n",
    "        else:\n",
    "            if count > max_count:\n",
    "                max_count = count\n",
    "                idx = i\n",
    "            count = 0\n",
    "    return max_count, idx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5, 55)"
      ]
     },
     "execution_count": 160,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max_count, pos = find_runs(tosses)\n",
    "max_count, pos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['T', 'H', 'H', 'H', 'H', 'H', 'T'],\n",
       "      dtype='<U1')"
      ]
     },
     "execution_count": 161,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tosses[(pos-max_count-1):(pos+1)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2**. Simulate a series of 10 sets of tosses of a fair coin in which you only stop when you find 5 consecutive heads. Find the average number of coins tossed in the 10 sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import itertools as it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def expt(n=5):\n",
    "    count = 0\n",
    "    for i in it.count(1):\n",
    "        toss = np.random.choice(list('HT'))\n",
    "        if toss == 'H':\n",
    "            count += 1\n",
    "            if count == n:\n",
    "                return i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.3000000000000007"
      ]
     },
     "execution_count": 171,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean([expt() for i in range(10)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Direct method using negative binomial sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.800000000000001"
      ]
     },
     "execution_count": 172,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(5 + np.random.negative_binomial(5, 0.5, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3**. A doubly stochastic matrix is one where each row and each column sums up to 1, and each cell value can be interpreted as a probability (i.e. it is between 0 and 1). Construct a doubly stochastic matrix of size (4,4) starting from a (4,4) matrix of standard random uniform samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.random.random((4,4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 191,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rsum = None\n",
    "csum = None\n",
    "\n",
    "while (np.any(rsum != 1)) | (np.any(csum != 1)):\n",
    "    x /= x.sum(0)\n",
    "    x = x / x.sum(1)[:, np.newaxis]\n",
    "    rsum = x.sum(1)\n",
    "    csum = x.sum(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.04242593,  0.09297455,  0.597877  ,  0.26672252],\n",
       "       [ 0.51742067,  0.07877105,  0.19201241,  0.21179587],\n",
       "       [ 0.07656644,  0.50788136,  0.10556133,  0.30999087],\n",
       "       [ 0.36358697,  0.32037303,  0.10454926,  0.21149075]])"
      ]
     },
     "execution_count": 192,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.,  1.,  1.,  1.])"
      ]
     },
     "execution_count": 193,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.sum(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.,  1.,  1.,  1.])"
      ]
     },
     "execution_count": 194,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.sum(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}