Monte_Carlo_SP21.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Monte Carlo methods\n",
    "\n",
    "This notebook will guide you through a few examples of how the power and speed of the computer processor can shed light on certain statistics related to a wide range of problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup\n",
    "import numpy as np\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "from pandas import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1: Finding area of a region\n",
    "\n",
    "By randomly selecting points and counting how many are in the region and how many are outside the region, you can get a good estimate of the true area of the region."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will try to estimate the area of the unit circle. It is all of the points (x,y) such that $x^2+y^2\\leq 1$. We will generate random points in the square with corners (1,1),(1,-1),(-1,1),(-1,-1). If it falls within the circle, we will increment ```circle_points```. Then divide ```circle_points``` by ```N=1000```, the total number of points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "circle_points = 0\n",
    "for trial in range(N):\n",
    "    r_x = 2*random.random()-1\n",
    "    r_y = 2*random.random()-1\n",
    "    if r_x**2 + r_y**2 <= 1:\n",
    "        circle_points += 1\n",
    "print(circle_points/N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tells us that $\\frac{\\text{Area of circle}}{\\text{Area of square}}\\approx B$ where $B$ is the value you get from the experiment. Since we know that the area of the square is $2^2= 4$, the area of the circle must be $4*B$. Does your answer line up with the geometric definition of the area of a circle?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2: Findmax\n",
    "\n",
    "As we have seen in class, the algorithm to find the maximum of a list of $n$ distinct integers is expected to update the maximum $H(n)-1$ times where $H(n)$ is the $n$th Harmonic number, i.e. $H(n) = 1+1/2+1/3+\\dots + 1/n$. We will conduct the following experiment: generate a random permutation $q$ of length $n$, then run the findmax algorithm on $q$ and count the number of times the maximum was updated. We will conduct this experiment 1000 times and take the average. Then we will check that against the theoretical expected value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generating random permutations\n",
    "\n",
    "In order to generate random permutations, I will use the following unranking algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_permutation(n):\n",
    "    output = []\n",
    "    sorted_list = [i for i in range(n)]\n",
    "    for j in range(n):\n",
    "        r = random.randint(0,n-j-1)\n",
    "        next_element = sorted_list[r]\n",
    "        output = output + [next_element]\n",
    "        sorted_list.remove(next_element)\n",
    "    return output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use the code below to generate a random permutation of length 10 (with integers 0,1,...,9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RP = random_permutation(10)\n",
    "print(RP)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Findmax function\n",
    "\n",
    "Here we will write the findmax function that takes a list of $n\\geq 1$ integers and returns ```(max_value, num_max_update)``` where ```max_value``` is the value of the maximum and ```num_max_update``` is the number of times the maximum updates "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def findmax(n_list):\n",
    "    n = len(n_list)\n",
    "    num_max_update = 0\n",
    "    max_value = n_list[0]\n",
    "    for i in range(1,n):\n",
    "        if n_list[i] > max_value:\n",
    "            max_value = n_list[i]\n",
    "            num_max_update += 1\n",
    "    return(max_value,num_max_update)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try it out on some randomly generated inputs. (you can change the value of ```n``` to test other list lengths.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10\n",
    "RP = random_permutation(n)\n",
    "print(RP)\n",
    "print(findmax(RP))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Experiment\n",
    "\n",
    "Let's run this experiment ```N=10000``` times where ```n=10```. We will total up the number of times max is updated and divide out by ```N``` to estimate the expected value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10000\n",
    "n = 10\n",
    "total = 0\n",
    "for j in range(N):\n",
    "    RP = random_permutation(10)\n",
    "    FM = findmax(RP)\n",
    "    total = total + FM[1]\n",
    "E_exp = total/N\n",
    "print(E_exp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Harmonic numbers\n",
    "\n",
    "Let's compare our experiment with the theoretical expected value.\n",
    "Recall that the theoretical expected value for inputs of length $n$ is $H(n)-1$ where $H(n)$ is the $n$th harmonic number.\n",
    "First let's generate a list of the harmonic numbers from 0 up to $100$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = [0]*101\n",
    "for i in range(1,101):\n",
    "    H[i] = H[i-1] + 1/i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check to see that ```H[10] = 2.928...```,  ```H[50] = 4.499...``` and ```H[100] = 5.187...```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(H[10])\n",
    "print(H[50])\n",
    "print(H[100])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the experimental error by taking the absolute value of the difference and dividing by the theoretical value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(n)\n",
    "print(E_exp)\n",
    "print(H[n]-1)\n",
    "exp_error = abs(E_exp-(H[n]-1))/(H[n]-1)\n",
    "print(exp_error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How well does the experiment estimate the theoretical expected value? By making ```N``` larger, the experiment should give better and better estimates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 3 Collecting a set of baseball cards.\n",
    "\n",
    "Let's suppose that in each pack of bubble gum you get one baseball card. There are ```n = 10``` different baseball cards and you want to collect them all. What is the expected number of packs of bubble gum you must buy before you collect all 10?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Experiment\n",
    "\n",
    "We will initialize an array of 10 0's. Then we will generate a number randomly from 0 to 9. We will change the corresponding 0 to a 1 until the array has all 1's.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def baseball_cards(n):\n",
    "    card_array = [0]*n\n",
    "    packs = 0\n",
    "    while 0 in card_array:\n",
    "        r = random.randint(0,n-1)\n",
    "        card_array[r] = 1\n",
    "        packs += 1\n",
    "    return packs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10\n",
    "baseball_cards(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's run the experiment ```N=1000``` times and take the average."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "n = 10\n",
    "total_packs = 0\n",
    "for trial in range(N):\n",
    "    total_packs = total_packs + baseball_cards(n)\n",
    "print(total_packs/N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try this experiment for different values of n. See if you can make a conjecture about the theoretical expected value for an arbitrary n. We will calculate this in class monday."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 4: Random Graphs\n",
    "\n",
    "There are a few different ways to define a random graph. We will be looking at the Erdős–Rényi model which is a randomly generated undirected graph.\n",
    "\n",
    "A $G(n,p)$ graph is constructed by starting with $n$ vertices and for every pair of vertices, there is a $p$ probability that there is an edge between the vertices and $1-p$ probability of no edge.\n",
    "\n",
    "We will be encoding the graph using an adjacency matrix.\n",
    "\n",
    "There are many fascinating probability questions that you can ask about random graphs. Let's try to experimentally approximate the answer to the questions:\n",
    "\n",
    "\n",
    "* What is the probability that the graph is connected?\n",
    "* What is the expected size of the connected component starting from vertex 0?\n",
    "* What is the expected number of connected components?\n",
    "* What is the expected maximum size of connected components?\n",
    "* What is the probability that the graph has an Eulerian Path?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generating $G(n,p)$\n",
    "\n",
    "This function generates a random graph given parameters ```n``` and ```p```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_G_n_p(n,p):\n",
    "    graph = [[0]*n for k in range(n)]     # this generates an n by n matrix of all 0's\n",
    "    for i in range(n-1):\n",
    "        for j in range(i+1,n):\n",
    "            r = random.random()\n",
    "            if r<p:                      # r is a random number chosen from the interval [0,1]\n",
    "                graph[i][j] = 1          # since the graph is undirected, you need to assign the same entry\n",
    "                graph[j][i] = 1          # for (i,j) and (j,i)\n",
    "    return graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try out the function for different values of ```n``` and ```p```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "generate_G_n_p(10,1/5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Graphsearch algorithm\n",
    "\n",
    "This algorithm is a version of graph search we have seen in class. It will use an array ```Status``` to keep track of the sets $X,F,U$. Then ```F``` will be used as a queue where the first element is the next to be considered and the new elements will be added to the end."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Graphsearch(G,s):                   # G is an n by n adjacency matrix of a graph and s (starting vertex) is an integer from 0 to n-1\n",
    "    n = len(G)\n",
    "    Status = ['U']*n\n",
    "    Status[s] = ['F']                   # initialize the Status array\n",
    "    F = [s]                             # F is a queue with only the vertex s\n",
    "    while len(F)>0:\n",
    "        w = F[0]\n",
    "        for i in range(n):\n",
    "            if G[w][i] == 1:\n",
    "                if Status[i] == 'U':\n",
    "                    Status[i] == 'F'\n",
    "                    F = F + [i]         # add vertex i to the queue F\n",
    "        Status [w] = 'X'\n",
    "        F = F[1:]                       # remove vertex w from the queue F\n",
    "    return Status"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a random 10 by 10 graph with p = 1/5 and run Graphsearch starting at 0. If ```Status``` is all ```'X'``` then the graph is connected. Otherwise, it is not. Is the graph_1 you generated connected or not?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graph_1 = generate_G_n_p(10,1/5)\n",
    "print(np.matrix(graph_1))\n",
    "Graphsearch(graph_1,0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Experiment 1\n",
    "\n",
    "Let's set ```n=10``` and ```p=1/4```. What is the probability the the graph is connected? We will run this experiment ```N=1000``` times and increment the total everytime the graph is connected. Then average the total."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "n = 10\n",
    "p = 1/4\n",
    "total_connected_graphs = 0\n",
    "for trial in range(N):\n",
    "    graphnp = generate_G_n_p(n,p)\n",
    "    GS = Graphsearch(graphnp,0)\n",
    "    if 'U' not in GS:\n",
    "        total_connected_graphs += 1\n",
    "print(total_connected_graphs/N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "n = 20\n",
    "p = 1/10\n",
    "total_isolated_vertices = 0\n",
    "for trial in range(N):\n",
    "    graphnp = generate_G_n_p(n,p)\n",
    "    for j in range(n):\n",
    "        if sum(graphnp[j]) == 1:\n",
    "            total_isolated_vertices += 1\n",
    "print(total_isolated_vertices/N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happens when you increase ```n```, does the result change?\n",
    "\n",
    "What happens when you change ```p```?, how does the result change?\n",
    "\n",
    "for ```n=10``` let's plot all results for ```p=1/x, 2/x, 3/x, ..., x/x``` for ```x=10```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Experiment 2\n",
    "\n",
    "I want to know in general how are ```n```, ```p```, and connectedness all related. For n=10 let's compute all results for p=1/x, 2/x, 3/x, ..., (x-1)/x for x=20. Then plot them using matplotlib. This might take a while to compute for large values of n."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = 20\n",
    "n = 10\n",
    "N = 1000\n",
    "data_points = [0]*x\n",
    "for i in range(1,x):\n",
    "    total_connected_graphs = 0\n",
    "    for trial in range(N):\n",
    "        graphnp = generate_G_n_p(n,i/(x))\n",
    "        GS = Graphsearch(graphnp,0)\n",
    "        if 'U' not in GS:\n",
    "            total_connected_graphs += 1\n",
    "    data_points[i] = total_connected_graphs/N\n",
    "print(data_points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot([j/x for j in range(x)],data_points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will see a drastic jump from very low probability to very high probability. This is a phenomenon known as a sharp threshold and when n is large, the jump gets more and more steep. This means that for very large graphs, if p is above the threshhold, it is almost certainly connected and if p is below the threshold then it is almost certainly not connected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Experiment 3:\n",
    "what is the expected size of the component that 0 belongs in? We will set ```n=10``` and ```p=1/n=0.1```. To find the size of the component that 0 belongs in, just count the number of ```'X'```'s in the output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "n = 10\n",
    "p = 1/n\n",
    "component_size_total = 0\n",
    "for trial in range(N):\n",
    "    graphnp = generate_G_n_p(n,p)\n",
    "    GS = Graphsearch(graphnp,0)\n",
    "    component_size = GS.count('X')\n",
    "    component_size_total = component_size_total + component_size\n",
    "print(component_size_total/N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to find a relationship between the size of the component containing 0 and G(n,1/n), we can run the experiment for 5,10,15,20,...,100 and plot the results to see if there is a recognizable curve. (This will take several seconds.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 5000\n",
    "data_points = [0]*11\n",
    "for j in range(1,11):\n",
    "    n = j*10\n",
    "    p = 1/n\n",
    "    component_size_total = 0\n",
    "    for trial in range(N):\n",
    "        graphnp = generate_G_n_p(n,p)\n",
    "        GS = Graphsearch(graphnp,0)\n",
    "        component_size = GS.count('X')\n",
    "        component_size_total = component_size_total + component_size\n",
    "    data_points[j] = component_size_total/N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot([j*10 for j in range(11)],data_points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The graph below graphs the same data on a log-log scale. The slope of this graph tells you the exponent of n for the growth rate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### %matplotlib qt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot([np.log(j*10) for j in range(1,11)],[np.log(k) for k in data_points[1:]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 5: Colonel Blotto\n",
    "\n",
    "Suppose I want to know about knight arrangements in colonel blotto. Particularly, I want to know what is the expected maximum number of knights in a single castle. We can randomly choose a number between 0 and ${109\\choose 9}-1$ then unrank it into a binary string of length 109 with exactly 9 ones. Then convert it to a knight arrangement and calculate the maximum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a function to calculate the binomial coefficient:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def binomial(n,k):\n",
    "    if k<0:\n",
    "        return 0\n",
    "    pr = 1\n",
    "    for i in range(k):\n",
    "        pr = pr*(n-i)/(i+1)\n",
    "    return int(pr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a recursive form of the unranking algorithm we learned in class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def unrank_binomial(n,k,d):\n",
    "    if n == k:\n",
    "        return [1]*n\n",
    "    if k == 0:\n",
    "        return [0]*n\n",
    "    if d < binomial(n-1,k):\n",
    "        return [0] + unrank_binomial(n-1,k,d)\n",
    "    else:\n",
    "        return [1] + unrank_binomial(n-1,k-1,d-binomial(n-1,k))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a procedure that converts a fixed density binary string into a knight arrangement sequence of 10 castles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bs2castles(binarystring):\n",
    "    L = len(binarystring)\n",
    "    t = sum(binarystring)\n",
    "    output = [0]*(t+1)\n",
    "    output_index = 0\n",
    "    for i in binarystring:\n",
    "        if i == 0:\n",
    "            output[output_index] += 1\n",
    "        else:\n",
    "            output_index += 1\n",
    "    return output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try out the sequence of commands to generate a random knight arrangement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "upper_rank_bound = binomial(109,9)\n",
    "RI = random.randint(0,upper_rank_bound-1)\n",
    "print(RI)\n",
    "UR = unrank_binomial(109,9,RI)\n",
    "print(UR)\n",
    "Knight_arrangement = bs2castles(UR)\n",
    "print(Knight_arrangement)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Experiment\n",
    "We will generate ```N = 1000``` knight arrangement and average the total_max_knights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "total_max_knights = 0\n",
    "upper_rank_bound = binomial(109,9)\n",
    "for trial in range(1000):\n",
    "    RI = random.randint(0,upper_rank_bound-1)\n",
    "    UR = unrank_binomial(109,9,RI)\n",
    "    Knight_arrangement = bs2castles(UR)\n",
    "    max_knights = max(Knight_arrangement)\n",
    "    total_max_knights = total_max_knights + max_knights\n",
    "print(total_max_knights/N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 6: HW7 Q2\n",
    "\n",
    "In HW7 Q2, there is a game based on rolling a die. We pick a positive integer $r$ and roll $r$ times or until we roll a 1 (whichever comes first.) If we finish $r$ rolls then we add up all of the rolls to get our score. If we roll a 1 then we get a score of 0.\n",
    "\n",
    "We can model this by choosing a particular value for $r$ (let's start with $r=8$.) Then we will simulate a game by rolling the die over and over until we get to $8$ rolls or roll a 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "r = 8\n",
    "sum_of_scores = 0\n",
    "trial_count = 0\n",
    "for trial in range(N):\n",
    "    score = 0\n",
    "    end_turn = False\n",
    "    number_of_rolls = 1\n",
    "    while not end_turn:\n",
    "        RI = random.randint(1,6)\n",
    "        if RI > 1:\n",
    "            score = score + RI\n",
    "            number_of_rolls = number_of_rolls + 1\n",
    "        else:\n",
    "            score = 0\n",
    "            end_turn = True\n",
    "        if number_of_rolls > r:\n",
    "            end_turn = True\n",
    "    trial_count = trial_count + score\n",
    "print(trial_count/N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can make a list of all possible values of r from 1 to 10 and use that to get an idea of the optimal value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_points = [0]*10\n",
    "for r in range(10):\n",
    "    N = 100000\n",
    "    sum_of_scores = 0\n",
    "    trial_count = 0\n",
    "    for trial in range(N):\n",
    "        score = 0\n",
    "        end_turn = False\n",
    "        number_of_rolls = 1\n",
    "        while not end_turn:\n",
    "            RI = random.randint(1,6)\n",
    "            if RI > 1:\n",
    "                score = score + RI\n",
    "                number_of_rolls = number_of_rolls + 1\n",
    "            else:\n",
    "                score = 0\n",
    "                end_turn = True\n",
    "            if number_of_rolls > r:\n",
    "                end_turn = True\n",
    "        trial_count = trial_count + score\n",
    "    print(r,trial_count/N)\n",
    "    data_points[r] = trial_count/N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot([j for j in range(0,10)],data_points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if instead of choosing a particular number of rolls r until we stop, we chose a score target T such that we roll until our score is T or higher and then stop (unless we roll a 1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_points_T = [0]*30\n",
    "for T in range(30):\n",
    "    N = 100000\n",
    "    sum_of_scores = 0\n",
    "    trial_count = 0\n",
    "    for trial in range(N):\n",
    "        score = 0\n",
    "        end_turn = False\n",
    "        while not end_turn:\n",
    "            RI = random.randint(1,6)\n",
    "            if RI > 1:\n",
    "                score = score + RI\n",
    "            else:\n",
    "                score = 0\n",
    "                end_turn = True\n",
    "            if score >= T:\n",
    "                end_turn = True\n",
    "        trial_count = trial_count + score\n",
    "    print(T,trial_count/N)\n",
    "    data_points_T[T] = trial_count/N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot([j for j in range(0,30)],data_points_T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}