C3/w2/lab/C3_W2_Lab_2_Dice_Simulations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "56c31ec0",
   "metadata": {},
   "source": [
    "# Lab: Simulate Dice Throws with Numpy 🎲🤖\n",
    "\n",
    "Welcome! In this lab you will see how you can use numpy to simulate different scenarios involving rolling dices from rolling a single die up to doing multiple rolls and summing the result. You will also see how to handle situations in which one of the sides of the die is loaded (it has a greater probability of landing on that side comparing to the rest).\n",
    "\n",
    "Let's get started! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af4835c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0968af94",
   "metadata": {},
   "source": [
    "## Represent a die\n",
    "\n",
    "The first thing you will need is to define how many sides your dice will have. You can even go a step further and represent a die by using a numpy array and assigning to each side a label which will be equal to the number of that side:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70b01f08",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the desired number of sides (try changing this value!)\n",
    "n_sides = 6\n",
    "\n",
    "# Represent a die by using a numpy array\n",
    "die = np.array([i for i in range(1, n_sides+1)])\n",
    "\n",
    "die"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6be050b4",
   "metadata": {},
   "source": [
    "## Roll the die\n",
    "\n",
    "With your die ready it is time to roll it. For now you will assume that the die is fair, which means the probability of landing on each side is the same (it follows a uniform distribution). To achieve this behaviour you can use the function [np.random.choice](https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html), which given a numpy array returns one of the entries in it randomnly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "154f10d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run this cell multiple times (every time you should get a different result at random)\n",
    "np.random.choice(die)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12750521",
   "metadata": {},
   "source": [
    "This is great but if you wanted to roll the die 20 times you will need to run the cell 20 times and record each result. Now you need a way to simulate several rolls at the same time. For this you can define the number of rolls you desire and use a list comprehension to roll the die as many times as you like, you can also save every roll in a numpy array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6f516d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Roll the die 20 times\n",
    "n_rolls = 20\n",
    "\n",
    "# Save the result of each roll\n",
    "rolls = np.array([np.random.choice(die) for _ in range(n_rolls)])\n",
    "\n",
    "rolls"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acc775fb",
   "metadata": {},
   "source": [
    "Now you have a convenient way of keeping track of the result of each roll, nice!\n",
    "\n",
    "What is you would like to know the mean and variance of this process. For this you can use numpy's functions [np.mean](https://numpy.org/doc/stable/reference/generated/numpy.mean.html) and [np.var](https://numpy.org/doc/stable/reference/generated/numpy.var.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb5fbb69",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute mean of 20 rolls\n",
    "m = np.mean(rolls)\n",
    "\n",
    "# Compute variance of 20 rolls\n",
    "v = np.var(rolls)\n",
    "\n",
    "print(f\"mean of rolls: {m:.2f}\\nvariance of rolls: {v:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f8a9054",
   "metadata": {},
   "source": [
    "You can even check the distribution of the rolls by plotting a histogram of the numpy array that holds the result of each throw. For this you will use the plotting library Seaborn, concretely the [sns.histplot](https://seaborn.pydata.org/generated/seaborn.histplot.html) function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52fe0fd6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display histogram of 20 rolls\n",
    "sns.histplot(rolls)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d210eb7",
   "metadata": {},
   "source": [
    "You probably didn't get a distribution that looks uniform (since the results are random). This happened because you are only simulating 20 rolls so far. Now try doing the same but for 20000 rolls:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01b974f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_rolls = 20_000\n",
    "\n",
    "rolls = np.array([np.random.choice(die) for _ in range(n_rolls)])\n",
    "\n",
    "print(f\"mean of rolls: {np.mean(rolls):.2f}\\nvariance of rolls: {np.var(rolls):.2f}\")\n",
    "\n",
    "sns.histplot(rolls)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00c6a396",
   "metadata": {},
   "source": [
    "Does this plot and the metrics of mean and variance allign with what you have learned about the uniform distribution during the course?\n",
    "\n",
    "Simulations are a great way of contrasting results against analytical solutions. For example, in this case the theoretical mean and variance are 3.5 and 2.916 respectively (you can check the formulas to get this results [here](https://en.wikipedia.org/wiki/Discrete_uniform_distribution)). The important thing to keep in mind is that the more simulations you perform the closer your results will be to the analytical values so always choose an appropriate number of simulations! \n",
    "\n",
    "Numpy is quite fast so performing 20 thousand runs is done fairly quick."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "654b7935",
   "metadata": {},
   "source": [
    "## Summing the result of rolling twice\n",
    "\n",
    "Now you want to throw the die twice and record the sum of the two rolls. For this you can do as before and save all results of the first roll in a numpy array but this time you will have a second array that saves the results for the second rolls. \n",
    "\n",
    "To get the sum you can simply sum the two arrays. This is possible because numpy allows for vectorized operations such as this one. When you sum two numpy arrays you will get a new array that includes the element-wise sum of the elements in the arrays you summed up.\n",
    "\n",
    "Notice that now you can compute the the mean and variance for the first rolls, the second rolls and the sum of rolls. You can also compute the covariance between the first and second rolls:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dcce70f",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_rolls = 20_000\n",
    "\n",
    "# First roll (same as before)\n",
    "first_rolls = np.array([np.random.choice(die) for _ in range(n_rolls)])\n",
    "\n",
    "# Second roll (code is the same but saved in a new numpy array)\n",
    "second_rolls = np.array([np.random.choice(die) for _ in range(n_rolls)])\n",
    "\n",
    "# Sum both rolls (this is easy since numpy allows vectorization)\n",
    "sum_of_rolls = first_rolls + second_rolls\n",
    "\n",
    "# Print mean, variance and covariance\n",
    "print(f\"mean of first_rolls: {np.mean(first_rolls):.2f}\\nvariance of first_rolls: {np.var(first_rolls):.2f}\\n\")\n",
    "print(f\"mean of second_rolls: {np.mean(second_rolls):.2f}\\nvariance of second_rolls: {np.var(second_rolls):.2f}\\n\")\n",
    "print(f\"mean of sum_of_rolls: {np.mean(sum_of_rolls):.2f}\\nvariance of sum_of_rolls: {np.var(sum_of_rolls):.2f}\\n\")\n",
    "print(f\"covariance between first and second roll:\\n{np.cov(first_rolls, second_rolls)}\")\n",
    "\n",
    "# Plot histogram\n",
    "sns.histplot(sum_of_rolls, stat = \"probability\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bf4dbd4",
   "metadata": {},
   "source": [
    "The resulting plot looks pretty Gaussian, as you might expect. Notice that the covariance between the first and second rolls is very close to zero since these two processes are independant of one another.\n",
    "\n",
    "Also notice that you can change the stat displayed in the histogram by changing the `stat` parameter of the `sns.histplot` function. In the previous exercises you were displaying the frequency but in this latter one you are plotting the probability, which makes more sense in this context. To check what other stats are available you can check the [docs](https://seaborn.pydata.org/generated/seaborn.histplot.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f70efe53",
   "metadata": {},
   "source": [
    "## Using loaded dice\n",
    "\n",
    "So far you have only simulated dice that are fair (all of the sides on them have the same probability of showing up), but what about simulating loaded dice (one or more of the sides have a greater probability of showing up)?\n",
    "\n",
    "It is actually pretty simple. [np.random.choice](https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html) has support for these kind of scenarios by having a parameter `p` you can set. This parameter controls the probability of selecting each one of the entries in the array.\n",
    "\n",
    "To see it in action, code a function that returns the probabilities of the die landing on each side given that one of the sides must have twice as much probability as the rest of them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84478b71",
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_dice(n_sides, loaded_number):\n",
    "    \n",
    "    # All probabilities are initially the same\n",
    "    probs = np.array([1/(n_sides+1) for _ in range(n_sides)])\n",
    "    \n",
    "    # Assign the loaded side a probability that is twice as the other ones\n",
    "    probs[loaded_number-1] = 1 - sum(probs[:-1])\n",
    "    \n",
    "    # Check that all probabilities sum up to 1\n",
    "    if not np.isclose(sum(probs), 1):\n",
    "        print(\"All probabilities should add up to 1\")\n",
    "        return\n",
    "    \n",
    "    return probs "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45768e32",
   "metadata": {},
   "source": [
    "Before using this function, check how the probabilities of a fair die would look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c63f659b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute same probabilities for every side\n",
    "probs_fair_die = np.array([1/n_sides]*n_sides)\n",
    "\n",
    "# Plot probabilities\n",
    "fair_die_sides = sns.barplot(x=die, y=probs_fair_die)\n",
    "fair_die_sides.set_ylim(0,0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bbf1913",
   "metadata": {},
   "source": [
    "Now get the probabilities by using the `load_dice` function. Try changing the loaded side!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67536017",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get probabilities if die is loaded towards side 2\n",
    "probs_loaded_die = load_dice(n_sides, loaded_number=2)\n",
    "\n",
    "# Plot probabilities\n",
    "loaded_die_sides = sns.barplot(x=die, y=probs_loaded_die)\n",
    "loaded_die_sides.set_ylim(0,0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cdf0ec5",
   "metadata": {},
   "source": [
    "Now, feed the `probs_loaded_die` array into `np.random.choice` and see how this affect the metrics and plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d98459f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_rolls = 20_000\n",
    "\n",
    "# Notice that the p parameter is being set\n",
    "first_rolls = np.array([np.random.choice(die, p=probs_loaded_die) for _ in range(n_rolls)])\n",
    "\n",
    "second_rolls = np.array([np.random.choice(die, p=probs_loaded_die) for _ in range(n_rolls)])\n",
    "\n",
    "sum_of_rolls = first_rolls + second_rolls\n",
    "\n",
    "print(f\"mean of first_rolls: {np.mean(first_rolls):.2f}\\nvariance of first_rolls: {np.var(first_rolls):.2f}\\n\")\n",
    "print(f\"mean of second_rolls: {np.mean(second_rolls):.2f}\\nvariance of second_rolls: {np.var(second_rolls):.2f}\\n\")\n",
    "print(f\"mean of sum_of_rolls: {np.mean(sum_of_rolls):.2f}\\nvariance of sum_of_rolls: {np.var(sum_of_rolls):.2f}\\n\")\n",
    "print(f\"covariance between first and second roll:\\n{np.cov(first_rolls, second_rolls)}\")\n",
    "\n",
    "\n",
    "sns.histplot(sum_of_rolls, stat = \"probability\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5af6e690",
   "metadata": {},
   "source": [
    "Now the histogram is skewed towards some values since some sums are now more likely than others. Try changing the loaded side and see how the histogram changes!\n",
    "\n",
    "Notice that covariance is still very close to zero since there is not any dependance between rolls of the die."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06b40238",
   "metadata": {},
   "source": [
    "## Dependant Rolls\n",
    "\n",
    "To finish this lab you will now simulate the scenario in which the second roll depends on the result of the first one. Say that you are playing a variant of the game you have played so far and you only roll the die a second time if the result of the first roll is greater or equal to 4.\n",
    "\n",
    "Before doing the simulations reflect on what might happen in this scenario. Some behavior you will probably see:\n",
    "\n",
    "- 1 is now a possible result since if you get a 1 in the first roll you don't roll again\n",
    "- 1, 2 and 3 now have a greater chance of showing up\n",
    "- 4 is now not a possible result since you need to roll again if you get a 4 in the first roll\n",
    "\n",
    "To achieve this behaviour you can use the [np.where](https://numpy.org/doc/stable/reference/generated/numpy.where.html) function, which given a condition can be used to zero-out the elements that don't meet its criteria:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6ea0bed",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_rolls = 20_000\n",
    "\n",
    "first_rolls = np.array([np.random.choice(die) for _ in range(n_rolls)])\n",
    "\n",
    "second_rolls = np.array([np.random.choice(die) for _ in range(n_rolls)])\n",
    "\n",
    "# Preserve the result of the second throw only if the first roll was greater or equal to 4\n",
    "second_rolls = np.where(first_rolls>=4, second_rolls, 0)\n",
    "\n",
    "sum_of_rolls = first_rolls + second_rolls\n",
    "\n",
    "print(f\"mean of first_rolls: {np.mean(first_rolls):.2f}\\nvariance of first_rolls: {np.var(first_rolls):.2f}\\n\")\n",
    "print(f\"mean of second_rolls: {np.mean(second_rolls):.2f}\\nvariance of second_rolls: {np.var(second_rolls):.2f}\\n\")\n",
    "print(f\"mean of sum_of_rolls: {np.mean(sum_of_rolls):.2f}\\nvariance of sum_of_rolls: {np.var(sum_of_rolls):.2f}\\n\")\n",
    "print(f\"covariance between first and second roll:\\n{np.cov(first_rolls, second_rolls)}\")\n",
    "\n",
    "\n",
    "sns.histplot(sum_of_rolls, stat = \"probability\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "970303f1",
   "metadata": {},
   "source": [
    "Looks like all of the predictions of this new scenario indeed happened. Notice that the covariance now is nowhere near zero since there is a dependency between the first and the second roll!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b016cd1",
   "metadata": {},
   "source": [
    "**Now you have finished this ungraded lab, nice job!**"
   ]
  }
 ],
 "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}