{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " \n", "
\n", "\n", "# Exploratory Computing with Python\n", "*Developed by Mark Bakker*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Notebook 9: Discrete random variables\n", "In this Notebook you learn how to deal with discrete random variables. Many of the functions we will use are included in the `random` subpackage of `numpy`. We will import this package and call it `rnd` so that we don't have to type `np.random.` all the time." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import numpy.random as rnd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Random numbers\n", "A random number generator lets you draw, at random, a number from a specified distribution. Several random number generators are included in the `random` subpackage of `numpy`. For example, the `randint(low, high, size)` function returns an integer array of shape `size` at random from `low` up to (but not including) `high`. For example, let's flip a coin 10 times and assign a 0 to heads and a 1 to tails. Note that the `high` is specified as `1 + 1`, which means it is `1` higher than the value we want." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0, 1, 1, 0, 1, 1, 1, 0, 1, 0])" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rnd.randint(0, 1 + 1, 10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we call the `randint` function again, we get a different sequence of heads (zeros) and tails (ones):" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0, 0, 1, 0, 1, 0, 1, 1, 0, 1])" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rnd.randint(0, 1+1, 10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Internally, the random number generator starts with what is called a *seed*. The seed is a number and is generated automatically (and supposedly at random) when you call the random number generator. The value of the seed exactly defines the sequence of random numbers that you get (so some people may argue that the generated sequence is at best pseudo-random, and you may not want to use the sequence for any serious cryptographic use, but for our purposes they are random enough). For example, let's set `seed` equal to 10" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1, 1, 0, 1, 0, 1, 1, 0, 1, 1])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rnd.seed(10)\n", "rnd.randint(0, 1 + 1, 10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we now specify the seed again as 10, we can generate the exact same sequence" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1, 1, 0, 1, 0, 1, 1, 0, 1, 1])" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rnd.seed(10)\n", "rnd.randint(0, 1+1, 10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The ability to generate the exact same sequence is useful during code development. For example, by seeding the random number generator, you can compare your output to output of others trying to solve the same problem. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Flipping a coin\n", "Enough for now about random number generators. Let's flip a coin 100 times and count the number of heads (0-s) and the number of tails (1-s):" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "number of heads: 59\n", "number of tails: 41\n" ] } ], "source": [ "flip = rnd.randint(0, 1 + 1, 100)\n", "headcount = 0\n", "tailcount = 0\n", "for i in range(100):\n", " if flip[i] == 0:\n", " headcount += 1\n", " else:\n", " tailcount += 1\n", "print('number of heads:', headcount)\n", "print('number of tails:', tailcount)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First of all, note that the number of heads and the number of tails adds up to 100. Also, note how we counted the heads and tails. We created counters `headcount` and `tailcount`, looped through all flips, and added 1 to the appropriate counter. Instead of a loop, we could have used a condition for the indices combined with a summation as follows" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "headcount 59\n", "tailcount 41\n" ] } ], "source": [ "headcount = np.count_nonzero(flip == 0)\n", "tailcount = np.count_nonzero(flip == 1)\n", "print('headcount', headcount)\n", "print('tailcount', tailcount)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How does that work? You may recall that the `flip == 0` statement returns an array with length 100 (equal to the lenght of `flip`) with the value `True` when the condition is met, and `False` when the condition is not met. The boolean `True` has the value 1, and the boolean `False` has the value 0. So we simply need to count the nonzero values using the `np.count_nonzero` function to find out how many items are `True`. \n", "\n", "The code above is easy, but if we do an experiment with more than two outcomes, it may be cumbersome to count the non-zero items for every possible outcome. So let's try to rewrite this part of the code using a loop. For this specific case, the number of lines of code doesn't decrease, but when we have an experiment with many different outcomes this will be much more efficient. Note that `dtype='int'` sets the array to integers." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "outcome 0 is 59\n", "outcome 1 is 41\n" ] } ], "source": [ "outcomes = np.zeros(2, dtype='int') # Two outcomes. heads are stored in outcome[0], tails in outcome[1]\n", "for i in range (2):\n", " outcomes[i] = np.count_nonzero(flip == i)\n", " print(f'outcome {i} is {outcomes[i]}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 1. Throwing a dice\n", "Throw a dice 100 times and report how many times you throw 1, 2, 3, 4, 5, and 6. Use a seed of 33. Make sure that the reported values add up to 100. Make sure you use a loop in your code as we did in the previous code cell." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Answers to Exercise 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Flipping a coin twice\n", "Next, we are going to flip a coin twice for 100 times and count the number of tails in two throws. We generate a random array of 0-s (heads) and 1-s (tails) with two rows (representing two coin flips) and 100 colums. The sum of the two rows represents the number of tails. Recall that the `np.sum` function takes an array as input argument and by default sums all the values in the array and returns one number. In this case we want to sum the rows. For that, the `sum` function has a keyword argument called `axis`, where `axis=0` sums over index 0 of the array (the rows), `axis=1` sums over the index 1 of the array (the columns), etc." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "number of 0, 1, 2 tails: [27 47 26]\n" ] } ], "source": [ "rnd.seed(55)\n", "flips = rnd.randint(low=0, high=1 + 1, size=(2, 100))\n", "tails = np.sum(flips, axis=0)\n", "number_of_tails = np.zeros(3, dtype='int')\n", "for i in range(3):\n", " number_of_tails[i] = np.count_nonzero(tails == i)\n", "print('number of 0, 1, 2 tails:', number_of_tails)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another way to simulate flipping a coin twice, is to draw a number at random from a set of 2 numbers (0 and 1). You need to replace the number after every draw, of course. The `numpy` function to draw a random number from a given array is called `choice`. The `choice` function has a keyword to specify whether values are replaced or not. Hence the following two ways to generate 5 flips are identical." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rnd.seed(55)\n", "flips1 = rnd.randint(low=0, high=1 + 1, size=5)\n", "rnd.seed(55)\n", "flips2 = rnd.choice(range(2), size=5, replace=True)\n", "np.alltrue(flips1 == flips2) # Check whether all values in the two arrays are equal" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Bar graph\n", "The outcome of the experiment may also be plotted with a bar graph" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAASp0lEQVR4nO3de9AddX3H8fdHLqJchEB0GCEGOwwjdRBqBERrBW8o1CsolipabLTeYu0FbC2Klwr1Rlsrlgo1KAoWUUScKlUirRckARQkWiiCBBgCQrh1xCZ8+8fZyEN4cp5NyJ7zPM++XzNnzu6e3T3fx4Of88vv/Pa3qSokSf3xiHEXIEkaLYNfknrG4JeknjH4JalnDH5J6pnNx11AGzvttFPNnz9/3GVI0oyybNmy26pq7rrbZ0Twz58/n6VLl467DEmaUZJcP9l2u3okqWcMfknqGYNfknrG4JeknjH4JalnDH5J6hmDX5J6xuCXpJ4x+CWpZ2bElbvql/nHnj/uEmat6044ZNwlaBqwxS9JPWPwS1LPGPyS1DMGvyT1jMEvST1j8EtSzxj8ktQzBr8k9YzBL0k9Y/BLUs8Y/JLUMwa/JPWMwS9JPWPwS1LPGPyS1DMGvyT1jMEvST1j8EtSzxj8ktQzBr8k9YzBL0k9Y/BLUs8Y/JLUMwa/JPWMwS9JPWPwS1LPdB78STZLclmSrzXruyW5OMnVSc5KsmXXNUiSHjCKFv8iYPmE9ROBj1fV7sAdwNEjqEGS1Og0+JPsAhwCfLpZD3AQcHazy2LgpV3WIEl6sK5b/CcBfwnc36zvCKyqqtXN+grg8R3XIEmaoLPgT3IosLKqlk3cPMmutZ7jFyZZmmTprbfe2kmNktRHXbb4nwG8OMl1wJkMunhOArZPsnmzzy7ATZMdXFWnVNWCqlowd+7cDsuUpH7pLPir6l1VtUtVzQeOAL5dVUcCFwKHNbsdBZzbVQ2SpIcaxzj+Y4B3JrmGQZ//qWOoQZJ6a/Opd3n4qmoJsKRZvhbYdxTvK0l6KK/claSeMfglqWcMfknqGYNfknrG4JeknjH4JalnDH5J6pkpgz/J4Um2bZbfneScJL/TfWmSpC60afH/TVXdneSZwAsYTKV8crdlSZK60ib41zTPhwAnV9W5gHfNkqQZqk3w35jkn4FXAl9P8siWx0mSpqE2Af5K4BvAwVW1CpgD/EWnVUmSOrPeSdqSzJmwumTCtvuApd2WJUnqyrDZOZcxuDtWJjyvVcATO6xLktSR9QZ/Ve02ykIkSaPRaj7+JDsAuwNbrd1WVRd1VZQkqTtTBn+SNwCLGNwf93Jgf+D7DO6hK0maYdqM6lkEPA24vqoOBPYBbu20KklSZ9oE/6+q6lcASR5ZVT8F9ui2LElSV9r08a9Isj3wFeCCJHcAN3VbliSpK1MGf1W9rFl8b5ILgccA/95pVZKkzgy7gGu7qrprnQu5rmietwFu77QySVInhrX4Pw8cyuQXcnkBlyTNUMMu4Do0SYDfq6pfjLAmSVKHho7qqaoCvjyiWiRJI9BmOOcPkjyt80okSSPRZjjngcAbk1wP3EvTx19Ve3VamSSpE22C/4WdVyFJGpk2XT0fqKrrJz6AD3RdmCSpG22C/7cnriTZDHhqN+VIkrq23uBP8q4kdwN7JbmredwNrATOHVmFkqRNar3BX1UfqqptgQ9X1XbNY9uq2rGq3jXCGiVJm9CUXT2GvCTNLm36+CVJs4jBL0k9M3QcfzNXz77A4xlMzHYT8MNmKoehkmwFXAQ8snmfs6vqPUl2A84E5gCXAq+pql8/rL9CktTasFE9zweuBt4LvAg4BDgeuLp5bSr3AQdV1VOAvYGDk+wPnAh8vKp2B+4Ajn5Yf4EkaYMMa/H/PfDcqrpu4samxf514EnDTtz8q+CeZnWL5lEMbtL+B832xQy+WE7ewLolSRtpWB//5sCKSbbfyCDEp5RksySXMxj7fwHwP8Cqqlrd7LKCQTeSJGlEhrX4TwMuSXImcEOzbVfgCODUNievqjXA3s09e7/M5P9KmPT3giQLgYUA8+bNa/N2kqQWhl7ABRzJYDbOpwMHNMtHNq+1VlWrgCXA/sD2SdZ+4ezCem7cXlWnVNWCqlowd+7cDXk7SdIQQ0f1VNVVwFXNfXerqu5oe+Ikc4H/q6pVSR4FPJfBD7sXAocxGNlzFE7/IEkjNWxUz7wkZyZZCVwM/DDJymbb/Bbn3hm4MMmPgUuAC6rqa8AxwDuTXAPsSMtuI0nSpjGsxX8WcBKDrp018JuZOQ9n0Frff9iJq+rHwD6TbL+WwbUBkqQxGDaqZ6eqOmtt6MPgx9qqOpNBS12SNAMNa/EvS/JJBmPtJ47qOQq4rOvCJEndGBb8r2VwVe3xDMbah8EXwHnYLy9JM9Z6g7+ZP+dkvKpWkmaVjZqdM8lxm7oQSdJobOy0zG/YpFVIkkZmvV09Se5a30vAo7opR5LUtWE/7q4CnlZVt6z7QpIbJtlfkjQDDAv+04EnAA8JfuDz3ZSz6c0/9vxxlzBrXXfCIeMuQdJGGDaq591DXjumm3IkSV3znruS1DMGvyT1zNBpmSVpKv6O1p2ufkdrFfzNrJyPm7h/Vf2ik4okSZ2aMviTvA14D4PRPfc3mwvYq8O6JEkdadPiXwTsUVW/7LoYSVL32vy4ewNwZ9eFSJJGo02L/1pgSZLzgfvWbqyqj3VWlSSpM22C/xfNY8vmIUmawaYM/qo6fhSFSJJGY9jsnCdV1TuSnMdgFM+DVNWLO61MktSJYS3+zzbPHxlFIZKk0Rg2Sduy5vk7oytHktQ15+qRpJ4x+CWpZ1oHf5KtuyxEkjQaUwZ/kgOSXAUsb9afkuSTnVcmSepEmxb/x4EXAL8EqKofAc/qsihJUndadfVU1bo3V1/TQS2SpBFoM2XDDUkOACrJlsDbabp9JEkzT5sW/5uAtwCPB1YAezfrkqQZqM1cPbcBR46gFknSCLQZ1bM4yfYT1ndIclq3ZUmSutKmq2evqlq1dqWq7gD26a4kSVKX2gT/I5LssHYlyRxa3qRdkjT9tAnwjwLfS3J2s3448MHuSpIkdWnKFn9VnQ4cBtwCrAReXlWfHX4UJNk1yYVJlif5SZJFzfY5SS5IcnXzvMNU55IkbTpt5+r5KXAOcC5wT5J5LY5ZDfxZVT0J2B94S5I9gWOBb1XV7sC3mnVJ0ohM2dWT5G3Aexi0+NcAYXBHrr2GHVdVNwM3N8t3J1nO4FqAlwDPbnZbDCwBjtmo6iVJG6xNH/8iYI+q+uXGvkmS+QxGAl0MPK75UqCqbk7y2PUcsxBYCDBvXpt/YEiS2mjT1XMDcOfGvkGSbYAvAe+oqrvaHldVp1TVgqpaMHfu3I19e0nSOtq0+K8FliQ5H7hv7caq+thUBybZgkHon1FV5zSbb0myc9Pa35nBD8aSpBFp0+L/BXABsCWw7YTHUEkCnAosX+dL4qvAUc3yUQx+MJYkjUibuXqOh8EduKrq3g049zOA1wBXJLm82fZXwAnAF5MczeBL5fANK1mS9HC0GdXzdAYt922AeUmeAryxqt487Liq+i8GI4Am85wNLVSStGm06eo5Ce/AJUmzhnfgkqSe8Q5cktQz3oFLknpmaIs/yWbAa6rKO3BJ0iwxtMVfVWsYzK0jSZol2vTxfzfJJ4CzgN+M46+qSzurSpLUmTbBf0Dz/L4J2wo4aNOXI0nqWpsrdw8cRSGSpNFoc+XucZNtr6r3TbZdkjS9tenqmTg/z1bAoTiOX5JmrDZdPR+duJ7kIwxm2JQkzUBt77k70aOBJ27qQiRJo9Gmj/8KBqN4ADYD5vLgET6SpBmkTR//oROWVwO3VNXqjuqRJHWsTVfPzsDtVXV9Vd0IbJVkv47rkiR1pE3wnwzcM2H9f5ttkqQZqE3wp6rW9vFTVffTrotIkjQNtQn+a5O8PckWzWMRcG3XhUmSutF2Pv4DgBsZzMe/H7Cwy6IkSd1pcwHXSuCIEdQiSRqBKVv8SRYn2X7C+g5JTuu2LElSV9p09exVVavWrlTVHcA+3ZUkSepSm+B/RJId1q4kmYOjeiRpxmoT4B8Fvpfk7Gb9cOCD3ZUkSepSmx93T0+ylAfuuPXyqrqq27IkSV1pOzvnFkAmLEuSZqg2o3oWAWcAOwGPBT6X5G1dFyZJ6kabPv6jgf2q6l6AJCcC3wf+scvCJEndaDVXD7BmwvoaHuj2kSTNMG1a/P8KXJzky836S4FTuytJktSlNqN6PpZkCfBMBi3911fVZV0XJknqRqsLsarqUuDSjmuRJI3AxtxsXZI0gxn8ktQznQV/ktOSrExy5YRtc5JckOTq5nmHYeeQJG16Xbb4PwMcvM62Y4FvVdXuwLeadUnSCHUW/FV1EXD7OptfAixulhczGBoqSRqhUffxP66qbgZonh+7vh2TLEyyNMnSW2+9dWQFStJsN21/3K2qU6pqQVUtmDt37rjLkaRZY9TBf0uSnQGa55Ujfn9J6r1RB/9XgaOa5aOAc0f8/pLUe10O5/wCg1k890iyIsnRwAnA85JcDTyvWZckjVBn986tqlev56XndPWekqSpTdsfdyVJ3TD4JalnDH5J6hmDX5J6xuCXpJ4x+CWpZwx+SeoZg1+Sesbgl6SeMfglqWcMfknqGYNfknrG4JeknjH4JalnDH5J6hmDX5J6xuCXpJ4x+CWpZwx+SeoZg1+Sesbgl6SeMfglqWcMfknqGYNfknrG4JeknjH4JalnDH5J6hmDX5J6xuCXpJ4x+CWpZwx+SeoZg1+Sesbgl6SeMfglqWcMfknqmbEEf5KDk/wsyTVJjh1HDZLUVyMP/iSbAf8EvBDYE3h1kj1HXYck9dU4Wvz7AtdU1bVV9WvgTOAlY6hDknpp8zG85+OBGyasrwD2W3enJAuBhc3qPUl+NoLapoOdgNvGXUQbOXHcFUwLM+bzAj+zxoz5zDbB5/WEyTaOI/gzybZ6yIaqU4BTui9nekmytKoWjLsOtePnNfP4mY2nq2cFsOuE9V2Am8ZQhyT10jiC/xJg9yS7JdkSOAL46hjqkKReGnlXT1WtTvJW4BvAZsBpVfWTUdcxjfWue2uG8/OaeXr/maXqId3rkqRZzCt3JalnDH5J6hmDf5pwGouZJclpSVYmuXLctWhqSXZNcmGS5Ul+kmTRuGsaJ/v4p4FmGov/Bp7HYLjrJcCrq+qqsRam9UryLOAe4PSqevK469FwSXYGdq6qS5NsCywDXtrX/4/Z4p8enMZihqmqi4Dbx12H2qmqm6vq0mb5bmA5g1kEesngnx4mm8ait/9RSl1KMh/YB7h4vJWMj8E/PbSaxkLSw5NkG+BLwDuq6q5x1zMuBv/04DQWUseSbMEg9M+oqnPGXc84GfzTg9NYSB1KEuBUYHlVfWzc9YybwT8NVNVqYO00FsuBLzqNxfSW5AvA94E9kqxIcvS4a9JQzwBeAxyU5PLm8aJxFzUuDueUpJ6xxS9JPWPwS1LPGPyS1DMGvyT1jMEvST1j8Kv3kixJ0vnNt5O8vZkd8ox1tu/dZmhhkgVJ/qFZfl2ST3RVq2a3kd96UZpNkmzeXIfRxpuBF1bVz9fZvjewAPj6sIOraimwdMOrlB7MFr9mhCTzm9byvzTzqX8zyaOa137TYk+yU5LrmuXXJflKkvOS/DzJW5O8M8llSX6QZM6Et/jDJN9LcmWSfZvjt27m3b+kOeYlE877b0nOA745Sa3vbM5zZZJ3NNs+BTwR+GqSP52w75bA+4BXNRcVvSrJvk0tlzXPezT7PjvJ1yZ5v8Ob9/pRkos2wf/cmuVs8Wsm2Z3BfQr+OMkXgVcAn5vimCczmIlxK+Aa4Jiq2ifJx4HXAic1+21dVQc08+yf1hz318C3q+qPkmwP/DDJfzT7Px3Yq6oeNDVzkqcCrwf2YzD53sVJvlNVb0pyMHBgVd22dv+q+nWS44AFVfXW5hzbAc+qqtVJngv8bfO3rs9xwAuq6samTmkog18zyc+r6vJmeRkwv8UxFzbzr9+d5E7gvGb7FcBeE/b7Agzm2U+yXROgzwdenOTPm322AuY1yxesG/qNZwJfrqp7AZKcA/wucFmbP7DxGGBxkt0ZzNK6xRT7fxf4TPNl2OvJx9SOXT2aSe6bsLyGBxouq3ngv+Wthhxz/4T1+3lww2fduUuKQYv9FVW1d/OYV1XLm9fvXU+Nk02xvaHez+AL68nA7/PQv+nBhVa9CXg3gxleL0+y4yaoQbOYwa/Z4Drgqc3yYRt5jlcBJHkmcGdV3clg0ry3NTM7kmSfFue5CHhpkkcn2Rp4GfCfUxxzN7DthPXHADc2y6+b6g2T/FZVXVxVxwG38eApvqWHMPg1G3wE+JMk3wN22shz3NEc/ylg7Uyb72fQzfLj5qbq75/qJM3t/T4D/JDBHZ4+XVVTdfNcCOy59sdd4O+ADyX5LrBZi9o/nOSKpsaLgB+1OEY95uycktQztvglqWcMfknqGYNfknrG4JeknjH4JalnDH5J6hmDX5J65v8BTCO52k6Jek8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.bar(range(0, 3), number_of_tails)\n", "plt.xticks(range(0, 3))\n", "plt.xlabel('number of tails')\n", "plt.ylabel('occurence in 100 trials');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Cumulative Probability\n", "Next, we compute the experimental probability of 0 tails, 1 tail, and 2 tails through division by the total number of trials (one trial is two coin flips). The three probabilities add up to 1. The cumulative probability distribution is obtained by cumulatively summing the probabilities using the `cumsum` function of `numpy`. The first value is the probability of throwing 0 tails. The second value is the probability of 1 or fewer tails, and the third value it the probability of 2 or fewer tails. The probability is computed as the number of tails divided by the total number of trials." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "cum_prob [0.27 0.74 1. ]\n" ] } ], "source": [ "prob = number_of_tails / 100 # number_of_tails was computed two code cells back\n", "cum_prob = np.cumsum(prob) # So cum_prob[0] = prob[0], cum_prob[1] = prob[0] + prob[1], etc.\n", "print('cum_prob ', cum_prob)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The cumulative probability distribution is plotted with a bar graph, making sure that all the bars touch each other (by setting the width to 1, in the case below)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAWEElEQVR4nO3debgldX3n8fcHEFFWpXtGZbFRUQfRAdNhTFwCLgngCESNQoIjBCXOBAEZ5xlmopiAE427RBKCS0SMIi7BFjBIEhbZtJudhhBbhKEDGVpFBIwg+J0/TjUcbp97bzV965y+Xe/X85zn1l7fe0/3+Zz6VdWvUlVIkvpro0kXIEmaLINAknrOIJCknjMIJKnnDAJJ6rlNJl3A2lqwYEEtWrRo0mVI0rxyxRVX/LCqFo6aN++CYNGiRSxbtmzSZUjSvJLk1unm2TQkST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs91FgRJPpPkziTXTzM/SU5MsiLJtUle2FUtkqTpdXlE8Flg7xnm7wPs3LwOB/6yw1okSdPoLAiq6iLgxzMssj/wuRq4HNgmyVO7qkeSNNok7yzeDrhtaHxlM+2OqQsmOZzBUQM77rjjWIqT+mzRsWdPugSNcMv7X93Jdid5sjgjpo18XFpVnVJVi6tq8cKFI7vKkCQ9RpMMgpXADkPj2wO3T6gWSeqtSQbBEuC/NFcPvQi4u6rWaBaSJHWrs3MESb4I7AksSLISeA/wOICqOhk4B9gXWAH8DDi0q1okSdPrLAiq6qBZ5hfwh13tX5LUjncWS1LPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk912kQJNk7yU1JViQ5dsT8HZOcn+SqJNcm2bfLeiRJa+osCJJsDJwE7APsAhyUZJcpi70LOKOqdgcOBP6iq3okSaN1eUSwB7Ciqm6uqgeA04H9pyxTwFbN8NbA7R3WI0kaocsg2A64bWh8ZTNt2B8DBydZCZwDvH3UhpIcnmRZkmWrVq3qolZJ6q0ugyAjptWU8YOAz1bV9sC+wGlJ1qipqk6pqsVVtXjhwoUdlCpJ/dVlEKwEdhga3541m34OA84AqKrLgM2ABR3WJEmaYtYgSPLVJK8e9U19FkuBnZPslGRTBieDl0xZ5v8Cr2j28x8YBIFtP5I0Rm0+3P8S+F3ge0nen+S5bTZcVQ8CRwDnAjcyuDpoeZLjk+zXLPbfgbcmuQb4InBIVU1tPpIkdWiT2Raoqr8H/j7J1gza9M9LchvwSeDzVfWLGdY9h8FJ4OFpxw0N3wC8+DHWLkmaA62ae5JsCxwCvAW4Cvg48ELgvM4qkySNxaxHBEm+BjwXOA14TVXd0cz6UpJlXRYnSererEEAfKpp4nlYksdX1f1VtbijuiRJY9Kmaei9I6ZdNteFSJImY9ojgiRPYXAn8BOS7M4jN4htBTxxDLVJksZgpqah32Jwgnh74CND0+8B/neHNUmSxmjaIKiqU4FTk7yuqr46xpokSWM0U9PQwVX1eWBRkmOmzq+qj4xYTVpri449e9IlSL02U9PQ5s3PLcZRiCRpMmZqGvqr5uefjK8cSdK4zdQ0dOJMK1bVkXNfjiRp3GZqGrpibFVIkiZmtquGJEkbuJmahj5WVUcn+QZrPlmMqtpvxGqSpHlmpqah05qfHxpHIZKkyZipaeiK5ueFzRPGnsvgyOCmqnpgTPVJkjrWphvqVwMnA99n0N/QTkn+oKq+2XVxkqTutemG+sPAXlW1AiDJM4GzAYNAkjYAbbqhvnN1CDRuBu7sqB5J0pjNdNXQa5vB5UnOAc5gcI7gd4ClY6hNkjQGMzUNvWZo+P8Bv9EMrwKe1FlFkqSxmumqoUPHWYgkaTLaXDW0GXAY8Dxgs9XTq+r3O6xLkjQmbU4WnwY8hcETyy5k8MSye7osSpI0Pm2C4FlV9W7gvqb/oVcDz++2LEnSuLQJgl80P3+SZFdga2BRZxVJksaqzQ1lpyR5EvBuYAmDJ5a9u9OqJEljM2sQVNWnmsELgWd0W44kadxmbRpKsm2SP09yZZIrknwsybbjKE6S1L025whOZ9ClxOuA1wM/BL7UZVGSpPFpc47gyVV1wtD4e5Mc0FVBkqTxanNEcH6SA5Ns1LzewKD3UUnSBmDaIEhyT5KfAn8AfAF4oHmdDryjzcaT7J3kpiQrkhw7zTJvSHJDkuVJvrD2v4IkaV3M1NfQluuy4SQbAycBrwJWAkuTLKmqG4aW2Rn4X8CLq+quJP9uXfYpSVp7bc4RkGQ/4GXN6AVVdVaL1fYAVlTVzc02Tgf2B24YWuatwElVdRdAVfmcA0kaszaXj74fOIrBB/gNwFHNtNlsB9w2NL6ymTbs2cCzk1yS5PIke7crW5I0V9ocEewL7FZVvwRIcipwFTCyzX9IRkyrEfvfGdiTQWd2306ya1X95FEbSg4HDgfYcccdW5QsSWqrzVVDANsMDW/dcp2VwA5D49sDt49Y5utV9Yuq+gFwE4NgeJSqOqWqFlfV4oULF7bcvSSpjTZB8D7gqiSfbY4GrgD+tMV6S4Gdk+yUZFPgQAZ9FQ07E9gLIMkCBk1FN7ctXpK07mZsGkoS4GLgRcCvMmju+Z9V9a+zbbiqHkxyBHAusDHwmapanuR4YFlVLWnm/WaSG4CHgP9RVT9ap99IkrRWZgyCqqokZ1bVr7Dmt/lZVdU5wDlTph03vH3gmOYlSZqANk1Dlyf51c4rkSRNRJurhvYC3pbkFuA+Bs1DVVUv6LIwSdJ4tAmCfTqvQpI0MW0eTHNrkhcCL2FwH8AlVXVl55VJksaizZ3FxwGnAtsCC4C/TvKurguTJI1Hm6ahg4Ddq+rn8HCXE1cC7+2yMEnSeLS5augWYLOh8ccD3++kGknS2LU5IrgfWJ7kPAbnCF4FXJzkRICqOrLD+iRJHWsTBH/bvFa7oJtSJEmT0OaqoVPHUYgkaTLa9j4qSdpAGQSS1HOtgyDJ5l0WIkmajDY3lP160030jc34f0zyF51XJkkaizZHBB8Ffgv4EUBVXcMjD7KXJM1zrZqGquq2KZMe6qAWSdIEtLmP4LYkvw5U88jJI2maiSRJ81+bI4K3AX8IbMfgYfO7NeOSpA1AmyOCVNXvdV6JJGki2hwRXJrkW0kOS7JN5xVJksZq1iCoqp2BdwHPA65MclaSgzuvTJI0Fm2vGvpuVR0D7AH8mMGDaiRJG4A2N5RtleTNSb4JXArcwSAQJEkbgDYni68BzgSOr6rLOq5HkjRmbYLgGVVVnVciSZqIaYMgyceq6mhgSZI1gqCq9uu0MknSWMx0RHBa8/ND4yhEkjQZ0wZBVV3RDO5WVR8fnpfkKODCLguTJI1Hm8tH3zxi2iFzXIckaUJmOkdwEPC7wE5JlgzN2pKmS2pJ0vw30zmC1fcMLAA+PDT9HuDaLouSJI3PTOcIbgVuBX5tfOVIksatzZ3FL0qyNMm9SR5I8lCSn46jOElS99qcLP4EcBDwPeAJwFuAP2+z8SR7J7kpyYokx86w3OuTVJLFbbYrSZo7bTudWwFsXFUPVdVfA3vNtk6SjYGTgH2AXYCDkuwyYrktGTz17DtrU7gkaW60CYKfNY+ovDrJB5K8A9i8xXp7ACuq6uaqegA4Hdh/xHInAB8Aft62aEnS3GkTBG8CNgaOAO4DdgBe12K97YDhh96vbKY9LMnuwA5VddZMG0pyeJJlSZatWrWqxa4lSW3N2ulcc/UQwL8Bf7IW286ozT08M9kI+Cgtbk6rqlOAUwAWL15sB3iSNIdmuqHsOoY+uKeqqhfMsu2VDI4eVtseuH1ofEtgV+CCJABPYdDB3X5VtWyWbUuS5shMRwT/eR23vRTYOclOwL8ABzK4UxmAqrqbwc1qACS5AHinISBJ4zXbDWWPWVU9mOQI4FwG5xg+U1XLkxwPLKuqJTNvQZI0DrOeI0hyD480EW0KPA64r6q2mm3dqjoHOGfKtOOmWXbP2bYnSZp7bU4Wbzk8nuQAfGaxJG0wWt1QNqyqzgRe3kEtkqQJaNM09Nqh0Y2AxcxwNZEkaX5p8/D61wwNPwjcwug7hCVJ81CbcwSHjqMQSdJktGka2gl4O7BoePmq2q+7siRJ49KmaehM4NPAN4BfdluOJGnc2gTBz6vqxM4rkSRNRJsg+HiS9wDfAu5fPbGqruysKknS2LQJgucz6Ir65TzSNFR4L4EkbRDaBMFvA89oHi4jSdrAtLmz+Bpgm64LkSRNRpsjgn8P/FOSpTz6HIGXj0rSBqBNELyn8yokSRPT5s7iC8dRiCRpMjp9HoEkaf3n8wgkqed8HoEk9ZzPI5CknvN5BJLUcz6PQJJ6rk3T0KnAUVX1k2b8ScCHq+r3uy5uri069uxJlyBJ6502J4tfsDoEAKrqLmD37kqSJI1TmyDYqDkKACDJk2l3bkGSNA+0+UD/MHBpkq8wuFroDcD/6bQqSdLYtDlZ/LkkyxjcOxDgtVV1Q+eVSZLGolUTT/PB74e/JG2A1vrOYknShsUgkKSeMwgkqecMAknqOYNAknqu0yBIsneSm5KsSHLsiPnHJLkhybVJ/iHJ07usR5K0ps6CIMnGwEnAPsAuwEFJdpmy2FXA4qp6AfAV4ANd1SNJGq3LI4I9gBVVdXNVPQCczpTuq6vq/Kr6WTN6ObB9h/VIkkboMgi2A24bGl/ZTJvOYcA3R81IcniSZUmWrVq1ag5LlCR1GQQZMW3kk82SHMzgyWcfHDW/qk6pqsVVtXjhwoVzWKIkqcteRFcCOwyNbw/cPnWhJK8E/gj4jaq6v8N6JEkjdHlEsBTYOclOSTYFDgSWDC+QZHfgr4D9qurODmuRJE2jsyCoqgeBI4BzgRuBM6pqeZLjk+zXLPZBYAvgy0muTrJkms1JkjrS6QNmquoc4Jwp044bGn5ll/uXJM3OO4slqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5zoNgiR7J7kpyYokx46Y//gkX2rmfyfJoi7rkSStqbMgSLIxcBKwD7ALcFCSXaYsdhhwV1U9C/go8Gdd1SNJGq3LI4I9gBVVdXNVPQCcDuw/ZZn9gVOb4a8Ar0iSDmuSJE2xSYfb3g64bWh8JfCfplumqh5McjewLfDD4YWSHA4c3ozem+SmTiqeXxYw5e+kifM9WT9tMO9L1q3N5OnTzegyCEZ9s6/HsAxVdQpwylwUtaFIsqyqFk+6Dj3C92T95Psyuy6bhlYCOwyNbw/cPt0ySTYBtgZ+3GFNkqQpugyCpcDOSXZKsilwILBkyjJLgDc3w68H/rGq1jgikCR1p7OmoabN/wjgXGBj4DNVtTzJ8cCyqloCfBo4LckKBkcCB3ZVzwbIprL1j+/J+sn3ZRbxC7gk9Zt3FktSzxkEktRzBsE8NFvXHRqvJJ9JcmeS6yddiwaS7JDk/CQ3Jlme5KhJ17Q+8xzBPNN03fHPwKsYXH67FDioqm6YaGE9luRlwL3A56pq10nXI0jyVOCpVXVlki2BK4AD/H8ymkcE80+brjs0RlV1Ed7/sl6pqjuq6spm+B7gRgY9GWgEg2D+GdV1h//ApWk0vRrvDnxnspWsvwyC+adVtxySIMkWwFeBo6vqp5OuZ31lEMw/bbrukHovyeMYhMDfVNXXJl3P+swgmH/adN0h9VrTnf2ngRur6iOTrmd9ZxDMM1X1ILC6644bgTOqavlkq+q3JF8ELgOek2RlksMmXZN4MfAm4OVJrm5e+066qPWVl49KUs95RCBJPWcQSFLPGQSS1HMGgST1nEEgST1nEGi9luSCJJ0/eDzJkU1PlX8zZfpubS47TLI4yYnN8CFJPtFyv09L8pW1rPXoJE9cm3VabPO5zSWWVyV5ZpJ7H2t9mn8MAm2wkqzNo1j/G7BvVf3elOm7AbMGQVUtq6oj16a+Zr3bq+r1a7na0cCcBgFwAPD1qtq9qr6/euJjrE/zjEGgdZZkUfNt+pNN3+/fSvKEZt7D3+iTLEhySzN8SJIzk3wjyQ+SHJHkmOYb6eVJnjy0i4OTXJrk+iR7NOtv3jwHYGmzzv5D2/1ykm8A3xpR6zHNdq5PcnQz7WTgGcCSJO8YWnZT4Hjgjc235Tcm2aOp5arm53OaZfdMctaI/f1Os69rklw0zd/u+qHav5bk75J8L8kHRix/JPA04Pymv/03JPlIM++oJDc3w89McnEz/Iqm3uuav9njp2xzXwbh8pYk589S39eb+m5K8p6h9+Ls5ne8Pskbp9at9VxV+fK1Ti9gEfAgsFszfgZwcDN8AbC4GV4A3NIMHwKsALYEFgJ3A29r5n2UQSdhq9f/ZDP8MuD6ZvhPh/axDYNnNGzebHcl8OQRdf4KcF2z3BbAcmD3Zt4twIIR6xwCfGJofCtgk2b4lcBXm+E9gbOmrtPsb7vVdU7zt7t+aL2bga2BzYBbgR1GrPNwrcBTgKXN8FcYdEGyHfBm4H3Ndm4Dnt0s87nVf9sp2/xj4J1D4/dOU98dwLbAE4DrgcXA61a/R81yW0/636SvtXt5RKC58oOquroZvoLBB8hszq+qe6pqFYMg+EYz/bop638RHu73f6sk2wC/CRyb5GoGYbEZsGOz/HlVNer5AC8B/raq7quqe4GvAS9t9+s9bGvgy8235I8Cz5tl+UuAzyZ5K7Bxi+3/Q1XdXVU/B24Anj7TwlX1r8AWzcNXdgC+wCAwXwp8G3gOg/fmn5tVTm3mP1bnVdWPqurfGPz9XsLg/Xplkj9L8tKqunsdtq8JMAg0V+4fGn4IWN0+/yCP/DvbbIZ1fjk0/suh9WHNbraLQXfcr6uq3ZrXjlV1YzP/vmlqHNWF99o6gUGA7Qq8hjV/p0cXWvU24F0MPqSvTrLtLNuf7u84k8uAQ4GbGHz4vxT4NQYhNBe/87A13osmZFYfbb0vyXFzvE91zCBQ125h8CEB8FhPOr4RIMlLgLubb5znAm9vepkkye4ttnMRcECSJybZHPhtBh+cM7mHQfPValsD/9IMHzLbDpM8s6q+U1XHAT/k0V2IP1ZTa7oIeGfz8ypgL+D+5u/0T8CiJM9qln0TcOE67PtVSZ7cnAM6ALgkydOAn1XV54EPAS9ch+1rAgwCde1DwH9NcimDcwSPxV3N+icDq3v2PAF4HHBt00xzwmwbqcGjCz8LfJfB06o+VVVXzbLa+cAuq08WAx9g8K33Eto19XywOUl7PYMP6mtarDObU4BvDp3Y/TaDgLmoqh5icE7gYoCmielQBs1Z1zE42jp5HfZ9MXAacDWD8yPLgOcD322a6f4IeO86bF8TYO+jklpJcgiDE/9HTLoWzS2PCCSp5zwikKSe84hAknrOIJCknjMIJKnnDAJJ6jmDQJJ67v8DhZdXR+bzXqIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.bar(range(0, 3), cum_prob, width=1)\n", "plt.xticks(range(0, 3))\n", "plt.xlabel('number of tails in two flips')\n", "plt.ylabel('cumulative probability');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 2. Flip a coin five times\n", "Flip a coin five times in a row and record how many times you obtain tails (varying from 0-5). Perform the experiment 1000 times. Make a bar graph with the total number of tails on the horizontal axis and the experimentally computed probability to get that many tails, on the vertical axis. Execute your code several times (hit [shift]-[enter]) and see that the graph changes a bit every time, as the sequence of random numbers changes every time. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute the cumulative probability. Print the values to the screen and make a plot of the cumulative probability function using a bar graph." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Answers to Exercise 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Probability of a Bernouilli variable\n", "In the previous exercise, we computed the probability of a certain number of heads in five flips experimentally. But we can, of course, compute the value exactly by using a few simple formulas. Consider the random variable $Y$, which is the outcome of an experiment with two possible values 0 and 1. Let $p$ be the probability of success, $p=P(Y=1)$. \n", "Then $Y$ is said to be a Bernoulli variable. The experiment is repeated $n$ times and we define $X$ as the number of successes in the experiment. The variable $X$ has a Binomial Distribution with parameters $n$ and $p$. The probability that $X$ takes value $k$ can be computed as (see for example [here](http://en.wikipedia.org/wiki/Binomial_distribution))\n", "\n", "$$P(X=k) = \\binom{n}{k}p^k(1-p)^{n-k}$$\n", "\n", "The term $\\binom{n}{k}$ may be computed with the `comb` function, which needs to be imported from the `scipy.special` package." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 3. Flip a coin 5 times revisited\n", "Go back to the experiment where we flip a coin five times in a row and record how many times we obtain tails.\n", "Compute the theoretical probability for 0, 1, 2, 3, 4, and 5 tails and compare your answer to the probability computed from 1000 trials, 10000 trials, and 100000 trials (use a loop for these three sets of trials). Do you approach the theoretical value with more trials?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Answers to Exercise 3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 4. Maximum value of two dice throws\n", "Throw a dice two times and record the maximum value of the two throws. Use the `np.max` function to compute the maximum value. Like the `np.sum` function, the `np.max` function takes an array as input argument and an optional keyword argument named `axis`. Perform the experiment 1000 times and compute the probability that the highest value is 1, 2, 3, 4, 5, or 6. Make a graph of the cumulative probability distribution function using a bar graph." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Answers to Exercise 4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 5. Maximum value of two dice throws revisited\n", "Refer back to Exercise 4.\n", "Compute the theoretical value of the probability of the highest dice when throwing the dice twice (the throws are labeled T1 and T2, respectively). There are 36 possible outcomes for this experiment. Let $M$ denote the random variable corresponding to this experiment (this means for instance that $M=3$ when your first throw is a 2, and the second throw is a 3). All outcomes of $M$ can easily be written down, as shown in the following Table: \n", "\n", "| T1$\\downarrow$ T2$\\to$ | 1 | 2 | 3 | 4 | 5 | 6 |\n", "|----|---|---|---|---|---|---|\n", "| 1 | 1 | 2 | 3 | 4 | 5 | 6 |\n", "| 2 | 2 | 2 | 3 | 4 | 5 | 6 |\n", "| 3 | 3 | 3 | 3 | 4 | 5 | 6 |\n", "| 4 | 4 | 4 | 4 | 4 | 5 | 6 |\n", "| 5 | 5 | 5 | 5 | 5 | 5 | 6 |\n", "| 6 | 6 | 6 | 6 | 6 | 6 | 6 |\n", "\n", "\n", "Use the 36 possible outcomes shown in the Table to compute the theoretical probability of $M$ being 1, 2, 3, 4, 5, or 6. Compare the theoretical outcome with the experimental outcome for 100, 1000, and 10000 dice throws." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Answers to Exercise 5" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generate random integers with non-equal probabilities\n", "So far, we have generated random numbers of which the probability of each outcome was the same (heads or tails, or the numbers on a dice, considering the throwing device was \"fair\"). What now if we want to generate outcomes that don't have the same probability? For example, consider the case that we have a bucket with 4 blue balls and 6 red balls. When you draw a ball at random, the probability of a blue ball is 0.4 and the probability of a red ball is 0.6. A sequence of drawing ten balls, with replacement, may be generated as follows" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "balls: [0 0 0 0 1 1 1 1 1 1]\n", "drawing 10 balls with replacement\n", "drawing: [1 1 1 1 1 0 1 1 1 0]\n", "blue balls: 2\n", "red balls: 8\n", "drawing 10 balls without replacement\n", "drawing: [0 0 1 1 0 1 0 1 1 1]\n", "blue balls: 4\n", "red balls: 6\n" ] } ], "source": [ "balls = np.zeros(10, dtype='int') # zero is blue\n", "balls[4:] = 1 # one is red\n", "print('balls:', balls)\n", "print('drawing 10 balls with replacement')\n", "rnd.seed(77)\n", "drawing = rnd.choice(balls, 10, replace=True)\n", "print('drawing:', drawing)\n", "print('blue balls:', np.count_nonzero(drawing == 0))\n", "print('red balls:', np.count_nonzero(drawing == 1))\n", "print('drawing 10 balls without replacement')\n", "rnd.seed(77)\n", "drawing = rnd.choice(balls, 10, replace=False)\n", "print('drawing:', drawing)\n", "print('blue balls:', np.count_nonzero(drawing == 0))\n", "print('red balls:', np.count_nonzero(drawing == 1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the previous example, we generated an example with 4 blue balls (zeros) and 6 red balls (ones), which was easy. But you can see this gets cumbersome when we have 4 million blue balls and 6 million red balls. But luckily, the `choice` function has an alternative: you can specifiy a sequence of values, and for each value you can specify the probability that it is drawn. Obviously, the `replace=False` argument doesn't make much sense anymore when the probabilities are specified (as the probabilities change when balls are not replaced). Repeating the previous example with replacement by specifying probabilities gives" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "drawing: [1 0 1 0 1 1 1 1 1 0]\n", "blue balls: 3\n", "red balls: 7\n" ] } ], "source": [ "drawing = rnd.choice([0, 1], 10, p=[0.4, 0.6]) # replace=False by default\n", "print('drawing:', drawing)\n", "print('blue balls:', np.count_nonzero(drawing == 0))\n", "print('red balls:', np.count_nonzero(drawing == 1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another cool feature of the `choice` function is that it also works on lists. For example, to pick 3 names at random from a list of 11 gives:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(['Daley', 'Quincy', 'Donny'], dtype='Election poll\n", "Consider an election where one million people will vote. 490,000 people will vote for candidate $A$ and 510,000 people will vote for candidate $B$. One day before the election, the company of 'Maurice the Dog' conducts a poll among 1000 randomly chosen voters. Compute whether the Dog will predict the winner correctly using the approach explained above and a seed of 2." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Perform the poll 1000 times. Count how many times the outcome of the poll is that candidate $A$ wins and how many times the outcome of the poll is that candidate $B$ wins. What is the probability that the Dog will predict the correct winner based on these 1000 polls of 1000 people? " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute the probability that the Dog will predict the correct winner based on 1000 polls of 5000 people. Does the probability that The Dog predicts the correct winner increase significantly when he polls 5000 people?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Answers to Exercise 6" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Answers to the exercises" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Answers to Exercise 1" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "number of times 1 is 17\n", "number of times 2 is 17\n", "number of times 3 is 15\n", "number of times 4 is 24\n", "number of times 5 is 19\n", "number of times 6 is 8\n", "total number of throws 100\n" ] } ], "source": [ "rnd.seed(33)\n", "dicethrow = rnd.randint(1, 6 + 1, 100)\n", "side = np.zeros(6, dtype='int')\n", "for i in range(6):\n", " side[i] = np.count_nonzero(dicethrow == i + 1)\n", " print(f'number of times {i + 1} is {side[i]}')\n", "print('total number of throws ', sum(side))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Back to Exercise 1\n", "\n", "Answers to Exercise 2" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAXOElEQVR4nO3df7BfdX3n8efLUMCiUpTsjgvEQM06DasFvQ3bgrS1iLi4YKcosOsKU3azVlm1jNPBsYMatlsUf3StdAU1BV1rivhjYk2LjIKKFkwQFAKLRIwS4wwoLGBV2MB7//ie2C9fPrn3e5N78r25PB8z3/me8znnc77vcyH3dc+P7+ekqpAkadSTJl2AJGl+MiAkSU0GhCSpyYCQJDUZEJKkpr0mXcBcOfDAA2vp0qWTLkOS9ig33HDDj6pqcWvZggmIpUuXsmHDhkmXIUl7lCTf29EyTzFJkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaFsw3qaU9zdJzPzfpEsay+YITJ12CJsQjCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUlOvAZHkhCS3J9mU5NzG8tckuTnJTUmuTbJ8aNmbu363J3lJn3VKkh6vt4BIsgi4CHgpsBw4fTgAOn9TVc+tqiOAdwLv6fouB04DDgdOAP6q254kaTfp8whiBbCpqu6sqoeBNcDJwytU1QNDs/sB1U2fDKypqoeq6rvApm57kqTdpM/B+g4C7hqa3wIcNbpSktcB5wB7Ay8a6nvdSN+DGn1XAisBlixZMidFS5IG+gyINNrqcQ1VFwEXJfkPwJ8CZ8yi7yXAJQBTU1OPW66FxdFPpd2rz1NMW4BDhuYPBrZOs/4a4OU72VeSNMf6DIj1wLIkhybZm8FF57XDKyRZNjR7InBHN70WOC3JPkkOBZYBX++xVknSiN5OMVXVtiRnA1cCi4DVVbUxySpgQ1WtBc5Ochzw/4D7GJxeolvvcuBWYBvwuqp6pK9aJUmP1+sT5apqHbBupO28oek3TNP3z4A/6686SdJ0/Ca1JKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWrqNSCSnJDk9iSbkpzbWH5OkluTfCvJF5I8a2jZI0lu6l5r+6xTkvR4e/W14SSLgIuAFwNbgPVJ1lbVrUOr3QhMVdVPk/wR8E7g1G7Zz6rqiL7qkyRNr88jiBXApqq6s6oeBtYAJw+vUFVXV9VPu9nrgIN7rEeSNAt9BsRBwF1D81u6th05C/j7ofl9k2xIcl2Sl/dRoCRpx3o7xQSk0VbNFZNXAVPAbw81L6mqrUkOA76Y5Oaq+s5Iv5XASoAlS5bMTdWSJKDfI4gtwCFD8wcDW0dXSnIc8BbgpKp6aHt7VW3t3u8ErgGOHO1bVZdU1VRVTS1evHhuq5ekJ7g+A2I9sCzJoUn2Bk4DHnM3UpIjgYsZhMPdQ+0HJNmnmz4QOBoYvrgtSepZb6eYqmpbkrOBK4FFwOqq2phkFbChqtYCFwJPAT6RBOD7VXUS8GvAxUkeZRBiF4zc/SRJ6lmf1yCoqnXAupG284amj9tBv68Bz+2zNknS9PwmtSSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqGisgkizquxBJ0vwy7hHEpiQXJlk+m40nOSHJ7Uk2JTm3sfycJLcm+VaSLyR51tCyM5Lc0b3OmM3nSpJ23bgB8Tzg28CHklyXZGWSp03XoTvquAh4KbAcOL0RMDcCU1X1POAK4J1d36cDbwWOAlYAb01ywJi1SpLmwFgBUVUPVtUHq+q3gD9h8Mv7h0kuS/LsHXRbAWyqqjur6mFgDXDyyHavrqqfdrPXAQd30y8Brqqqe6vqPuAq4IRZ7ZkkaZeMfQ0iyUlJPg38T+DdwGHAZ4F1O+h2EHDX0PyWrm1HzgL+fjZ9uyOZDUk23HPPPePsiiRpTHuNud4dwNXAhVX1taH2K5Icu4M+abRVc8XkVcAU8Nuz6VtVlwCXAExNTTW3LUnaOeNeg3h1VZ01HA5JjgaoqtfvoM8W4JCh+YOBraMrJTkOeAtwUlU9NJu+kqT+jBsQ72u0/eUMfdYDy5IcmmRv4DRg7fAKSY4ELmYQDncPLboSOD7JAd3F6eO7NknSbjLtKaYkvwn8FrA4yTlDi54GTPvdiKraluRsBr/YFwGrq2pjklXAhqpaC1wIPAX4RBKA71fVSVV1b5LzGYQMwKqquncn9k+StJNmugaxN4Nf4HsBTx1qfwA4ZaaNV9U6Ri5iV9V5Q9PHTdN3NbB6ps+QJPVj2oCoqi8BX0pyaVV9bzfVJEmaB2Y6xfQXVfVG4P1JWncRndRbZZKkiZrpFNNHu/d39V2IJGl+mekU0w3d+5d2TzmSpPliplNMN7ODL7cBdGMoSZIWoJlOMb1st1QhSZp3ZjrF5J1LkvQENe03qZNc270/mOSB0ffdU6IkaRJmOoI4pnt/6nTrSZIWnnFHcyXJ84FjGFy0vraqbuytKknSxI37PIjzgMuAZwAHApcm+dM+C5MkTda4RxCnA0dW1c8BklwAfAP4730VJkmarHGH+94M7Ds0vw/wnTmvRpI0b8z0Rbm/ZHDN4SFgY5KruvkXA9f2X54kaVJmOsW0oXu/Afj0UPs1vVQjSZo3ZrrN9bLdVYgkaX4Z6yJ1kmXAnwPLGboWUVWH9VSXJGnCxr1I/dfA/wK2Ab8LfIR/HgpckrQAjRsQT66qLwCpqu9V1duAF/VXliRp0sb9HsTPkzwJuCPJ2cAPgH/RX1mSpEkbNyDeCPwy8HrgfAZHD2f0VZR23dJzPzfpEsay+YITJ12CpB0YKyCqaj1AdxTx+qp6sNeqJEkTN+5dTFMMLlQ/tZu/H/jD7Y8klSTwyHWhGfci9WrgtVW1tKqWAq9jEBjTSnJCktuTbEpybmP5sUm+kWRbklNGlj2S5KbutXbMOiVJc2TcaxAPVtVXts9U1bVJpj3NlGQRcBGDYTm2AOuTrK2qW4dW+z5wJvCmxiZ+VlVHjFmfJGmOzTQW0/O7ya8nuRj4OIOxmE5l5uE2VgCbqurObltrgJOBXwREVW3ulj26E7VLkno00xHEu0fm3zo0XTP0PQi4a2h+C3DUmHUB7JtkA4Mv511QVZ8ZXSHJSmAlwJIlS2axaUnSTGYai+l3d2HbaW1yFv2XVNXWJIcBX0xyc1U9ZojxqroEuARgampqNtuWJM1g3CfK7Z/kPUk2dK93J9l/hm5bgEOG5g8Gto5bWFVt7d7vZHA668hx+0qSdt1s7mJ6EHhl93qAme9iWg8sS3Jokr2B04Cx7kZKckCSfbrpA4GjGbp2IUnq37h3Mf1qVf3B0Pzbk9w0XYeq2tYNy3ElsAhYXVUbk6wCNlTV2iS/weA5EwcA/z7J26vqcODXgIu7i9dPYnANwoCQpN1o3ID4WZJjqupagCRHAz+bqVNVrQPWjbSdNzS9nsGpp9F+XwOeO2ZtkqQejBsQrwE+MnTd4T4ci0mSFrQZA6Ibf+k5VfXrSZ4GUFUP9F6ZJGmiZrxIXVWPAmd30w8YDpL0xDDuXUxXJXlTkkOSPH37q9fKJEkTNe41iD9k8CW31460+0xqSVqgxg2I5QzC4RgGQfEV4AN9FSVJmrxxA+IyBl+Oe183f3rX9so+ipIkTd64AfGcqvr1ofmrk3yzj4IkSfPDuBepb0zyb7fPJDkK+Go/JUmS5oNxjyCOAl6d5Pvd/BLgtiQ3A1VVz+ulOknSxIwbECf0WoUkad4ZKyCq6nt9FyJJml/GvQYhSXqCMSAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqanXgEhyQpLbk2xKcm5j+bFJvpFkW5JTRpadkeSO7nVGn3VKkh6vt4BIsgi4CHgpg0eWnp5k+chq3wfOBP5mpO/TgbcyGGZ8BfDWJAf0Vask6fH6PIJYAWyqqjur6mFgDXDy8ApVtbmqvgU8OtL3JcBVVXVvVd0HXIVDjkvSbtVnQBwE3DU0v6Vrm7O+SVYm2ZBkwz333LPThUqSHq/PgEijreayb1VdUlVTVTW1ePHiWRUnSZpenwGxBThkaP5gYOtu6CtJmgN9BsR6YFmSQ5PsDZwGrB2z75XA8UkO6C5OH9+1SZJ2k94Coqq2AWcz+MV+G3B5VW1MsirJSQBJfiPJFuAVwMVJNnZ97wXOZxAy64FVXZskaTcZ65nUO6uq1gHrRtrOG5pez+D0UavvamB1n/VJknbMb1JLkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpKZeAyLJCUluT7IpybmN5fsk+dtu+fVJlnbtS5P8LMlN3esDfdYpSXq8vfracJJFwEXAi4EtwPoka6vq1qHVzgLuq6pnJzkNeAdwarfsO1V1RF/1SZKm1+cRxApgU1XdWVUPA2uAk0fWORm4rJu+Avi9JOmxJknSmPoMiIOAu4bmt3RtzXWqahtwP/CMbtmhSW5M8qUkL+yxTklSQ2+nmIDWkUCNuc4PgSVV9eMkLwA+k+TwqnrgMZ2TlcBKgCVLlsxByZKk7fo8gtgCHDI0fzCwdUfrJNkL2B+4t6oeqqofA1TVDcB3gH89+gFVdUlVTVXV1OLFi3vYBUl64uozINYDy5IcmmRv4DRg7cg6a4EzuulTgC9WVSVZ3F3kJslhwDLgzh5rlSSN6O0UU1VtS3I2cCWwCFhdVRuTrAI2VNVa4MPAR5NsAu5lECIAxwKrkmwDHgFeU1X39lWrJOnx+rwGQVWtA9aNtJ03NP1z4BWNfp8EPtlnbZKk6flNaklSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqSmXm9z3ZMsPfdzky5hLJsvOHHSJUhPGE/03wseQUiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaeg2IJCckuT3JpiTnNpbvk+Rvu+XXJ1k6tOzNXfvtSV7SZ52SpMfrLSCSLAIuAl4KLAdOT7J8ZLWzgPuq6tnAe4F3dH2XA6cBhwMnAH/VbU+StJv0eQSxAthUVXdW1cPAGuDkkXVOBi7rpq8Afi9JuvY1VfVQVX0X2NRtT5K0m+zV47YPAu4amt8CHLWjdapqW5L7gWd07deN9D1o9AOSrARWdrM/SXL73JQ+Zw4EfjSXG8w75nJrs7bQ9gcW3j4ttP2BhbdP821/nrWjBX0GRBptNeY64/Slqi4BLpl9abtHkg1VNTXpOubKQtsfWHj7tND2BxbePu1J+9PnKaYtwCFD8wcDW3e0TpK9gP2Be8fsK0nqUZ8BsR5YluTQJHszuOi8dmSdtcAZ3fQpwBerqrr207q7nA4FlgFf77FWSdKI3k4xddcUzgauBBYBq6tqY5JVwIaqWgt8GPhokk0MjhxO6/puTHI5cCuwDXhdVT3SV609mrenv3bSQtsfWHj7tND2BxbePu0x+5PBH+ySJD2W36SWJDUZEJKkJgOiBzMNMbKnSbI6yd1Jbpl0LXMhySFJrk5yW5KNSd4w6Zp2VZJ9k3w9yTe7fXr7pGuaC0kWJbkxyd9Nupa5kGRzkpuT3JRkw6TrmYnXIOZYNyTIt4EXM7hddz1welXdOtHCdkGSY4GfAB+pqn8z6Xp2VZJnAs+sqm8keSpwA/DyPfy/UYD9quonSX4JuBZ4Q1VdN0PXeS3JOcAU8LSqetmk69lVSTYDU1U1p1+U64tHEHNvnCFG9ihV9WUGd5ktCFX1w6r6Rjf9IHAbjW/q70lq4Cfd7C91rz36r78kBwMnAh+adC1PVAbE3GsNMbJH//JZyLoRhI8Erp9sJbuuOx1zE3A3cFVV7en79BfAnwCPTrqQOVTA55Pc0A0VNK8ZEHNvrGFCNHlJngJ8EnhjVT0w6Xp2VVU9UlVHMBh5YEWSPfZ0YJKXAXdX1Q2TrmWOHV1Vz2cwyvXrutO385YBMfccJmQP0J2n/yTwsar61KTrmUtV9X+BaxgMlb+nOho4qTtnvwZ4UZL/PdmSdl1Vbe3e7wY+zTwfpdqAmHvjDDGiCeou6H4YuK2q3jPpeuZCksVJfqWbfjJwHPB/JlvVzquqN1fVwVW1lMG/oS9W1asmXNYuSbJfd1MESfYDjgfm9Z2BBsQcq6ptwPYhRm4DLq+qjZOtatck+Tjwj8BzkmxJctaka9pFRwP/icFfpTd1r3836aJ20TOBq5N8i8EfKVdV1YK4NXQB+ZfAtUm+yWBsuc9V1T9MuKZpeZurJKnJIwhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZENojJbkmSe8Pfk/y+m7U14+NtB8xzq2xSaaSvK+bPjPJ+8f83H+V5IpZ1vrCbiTXm5IcNNv+02z3wm67FyZ5W5I3de2rkhw3F5+h+am3R45K81WSvbrvq4zjtcBLq+q7I+1HMBhldN10natqAzDrYZ27b9yeMstu/xF4V1X9dTc/2/478l+BxVX1UJK3bW+sqvPmaPuapzyCUG+SLO3++v5g9xfo57tv+T7mCCDJgd2QCtv/yv5Mks8m+W6Ss5Oc0z0T4LokTx/6iFcl+VqSW5Ks6Prv1z2/Yn3X5+Sh7X4iyWeBzzdqPafbzi1J3ti1fQA4DFib5I+H1t0bWAWc2v21fmqSFV0tN3bvz+nW/Z3WswySvKL7rG8m+fIOfna3DNX+qST/kOSOJO9srP+fgVcC5yX52Ej/65McPrTuNUlesKOf1ch21wL7AdcnOXVk2aVJTummNyd5RwbPpPh6kmePs5+a56rKl69eXsBSYBtwRDd/OfCqbvoaBuPiAxwIbO6mzwQ2AU8FFgP3A6/plr2XwcB62/t/sJs+Frilm/4fQ5/xKwyezbFft90twNMbdb4AuLlb7ynARuDIbtlm4MBGnzOB9w/NPw3Yq5s+DvhkN/07wN+N9uk+76Dtde7gZ3fLUL87gf2BfYHvAYc0+lwKnNLo/8fA27vpZwLfnu5n1djuT4am3wa8qfF5m4G3dNOvHtrnaffT1/x+eQShvn23qm7qpm9g8ItrJldX1YNVdQ+DgPhs137zSP+Pwy+eV/G0biyi44FzMxj2+hoGv1CXdOtfVVWt51ocA3y6qv6pBs9U+BTwwvF27xf2Bz7R/dX+XuDwGdb/KnBpkv8CLBpj+1+oqvur6ufArcCzZlHb5cAruulXAp/opqf7We2Mjw+9/2Y3Pdv91DziNQj17aGh6UeAJ3fT2/jnU5z7TtPn0aH5R3ns/7Oj48QUg+HW/6Cqbh9ekOQo4J92UGNriPbZOp9BsP1+Bs+YuGa6lavqNV1NJwI3JTmiqn48TZfRn+PY/3ar6gdJfpzkecCpDK4pwA5+VrugRqd3Yj81j3gEoUnZzODUDuz8xdRTAZIcA9xfVfczGCTxvyVJt+zIMbbzZeDlSX45g1E2fx/4ygx9HmRwGmy7/YEfdNNnzvSBSX61qq6vwYXeH/HYIeL7sIbBw3f2r6qbu7ad+VlN59Sh93/strm791NzyIDQpLwL+KMkX2NwDWJn3Nf1/wCwfYTZ8xk8bvNb3eme82faSA0eP3opgxE2rwc+VFU3ztDtamD59ovUwDuBP0/yVcY7lXJhBg+vv4VBQH1zjD674goGw2ZfPtQ265/VDPZJcj3wBgbXPWD376fmkKO5Stpl3V1oU1X1o0nXornjEYQkqckjCElSk0cQkqQmA0KS1GRASJKaDAhJUpMBIUlq+v9bxlH+MH/tZAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "N = 1000\n", "tails = np.sum(rnd.randint(0, 1 + 1, (5, 1000)), axis=0)\n", "counttails = np.zeros(6, dtype='int')\n", "for i in range(6):\n", " counttails[i] = np.count_nonzero(tails == i)\n", "plt.bar(range(0, 6), counttails / N)\n", "plt.xlabel('number of tails in five flips')\n", "plt.ylabel('probability');" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "cumprob: [0.033 0.192 0.491 0.812 0.968 1. ]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAW/UlEQVR4nO3de7hddX3n8ffHIGKRi5pMdbgYrLQOWgc0pbZaq1YtiIIXFGhpxUHRGVGUcWboVFHRaa13rYwWFUG0It6YKLFCW9B6AZNwk4RSU4xDxBniDQFHEPzOH3sFNyf7nLNyctbeOWe9X8+zn73u+7sOZH/2+q21fitVhSSpv+416QIkSZNlEEhSzxkEktRzBoEk9ZxBIEk9t9OkC9hWS5cureXLl0+6DElaUNauXfv9qlo2at6CC4Lly5ezZs2aSZchSQtKku9MN8+mIUnqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6rrMgSHJmkpuSXDPN/CR5T5INSa5O8uiuapEkTa/LI4KzgENmmH8osH/zOgF4X4e1SJKm0VkQVNWXgR/OsMgRwEdq4FJgzyQP7qoeSdJok7yzeC/ghqHxTc20701dMMkJDI4a2HfffcdSnKTFYfkpF0y6hHmz8c2HdbLdSZ4szohpIx+XVlVnVNWKqlqxbNnIrjIkSXM0ySDYBOwzNL43cOOEapGk3ppkEKwE/rS5euixwM1VtVWzkCSpW52dI0jyceCJwNIkm4DXAfcGqKr3A6uApwMbgJ8CL+yqFknS9DoLgqo6Zpb5Bbysq8+XJLWz4J5HIKl7i+lKG83OLiYkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6rtMgSHJIkuuSbEhyyoj5+ya5OMkVSa5O8vQu65EkbW2nrjacZAlwOvBUYBOwOsnKqlo/tNhrgPOq6n1JDgBWAcu7qknq0vJTLph0CdKcdHlEcDCwoaqur6o7gHOBI6YsU8DuzfAewI0d1iNJGqHLINgLuGFofFMzbdjrgWOTbGJwNPDyURtKckKSNUnWbN68uYtaJam3ugyCjJhWU8aPAc6qqr2BpwPnJNmqpqo6o6pWVNWKZcuWdVCqJPVXl0GwCdhnaHxvtm76OR44D6Cqvg7sAiztsCZJ0hSzBkGSTyc5bNQv9VmsBvZPsl+SnYGjgZVTlvnfwB80n/PvGASBbT+SNEZtvtzfB/wR8K0kb07y8DYbrqo7gROBLwLXMrg6aF2S05Ic3iz2n4EXJ7kK+DhwXFVNbT6SJHVo1stHq+rvgb9PsgeDNv2LktwAfAD4aFX9fIZ1VzE4CTw87dSh4fXA4+ZYuyRpHrRq7knyQOA44EXAFcC7gUcDF3VWmSRpLGY9IkjyGeDhwDnAM6vqe82sTyRZ02VxkqTutbmz+INNE8/dktynqm6vqhUd1SVJGpM2TUNvGjHt6/NdiCRpMqY9IkjyIAZ3At83yUH88gax3YFfGUNtkqQxmKlp6A8ZnCDeG3jH0PRbgP/eYU2SpDGaNgiq6mzg7CTPrapPj7EmSdIYzdQ0dGxVfRRYnuTkqfOr6h0jVpMkLTAzNQ3t2rzfbxyFSJImY6amob9p3t8wvnIkSeM2U9PQe2ZasapeMf/lSJLGbaamobVjq0KSNDGzXTUkSVrkZmoaeldVvTLJ59j6yWJU1eEjVpMkLTAzNQ2d07y/bRyFSJImY6amobXN+5eaJ4w9nMGRwXVVdceY6pMkdaxNN9SHAe8H/pVBf0P7JXlJVX2h6+IkSd1r0w3124EnVdUGgCS/BlwAGASStAi06Yb6pi0h0LgeuKmjeiRJYzbTVUPPaQbXJVkFnMfgHMHzgNVjqE2SNAYzNQ09c2j4/wK/3wxvBu7fWUWSpLGa6aqhF46zEEnSZLS5amgX4HjgEcAuW6ZX1X/osC5J0pi0OVl8DvAgBk8s+xKDJ5bd0mVRkqTxaRMED6uq1wK3Nf0PHQb8ZrdlSZLGpU0Q/Lx5/3GSRwJ7AMs7q0iSNFZtbig7I8n9gdcCKxk8sey1nVYlSRqbWYOgqj7YDH4JeGi35UiSxm3WpqEkD0zy10kuT7I2ybuSPHAcxUmSutfmHMG5DLqUeC5wJPB94BNdFiVJGp825wgeUFVvHBp/U5JndVWQJGm82hwRXJzk6CT3al7PZ9D7qCRpEZg2CJLckuQnwEuAvwXuaF7nAq9qs/EkhyS5LsmGJKdMs8zzk6xPsi7J3277LkiStsdMfQ3ttj0bTrIEOB14KrAJWJ1kZVWtH1pmf+DPgMdV1Y+S/Jvt+UxJ0rZrc46AJIcDT2hGL6mqz7dY7WBgQ1Vd32zjXOAIYP3QMi8GTq+qHwFUlc85kKQxa3P56JuBkxh8ga8HTmqmzWYv4Iah8U3NtGG/Dvx6kq8muTTJIe3KliTNlzZHBE8HDqyqXwAkORu4AhjZ5j8kI6bViM/fH3gig87s/inJI6vqx/fYUHICcALAvvvu26JkSVJbba4aAthzaHiPlutsAvYZGt8buHHEMv+rqn5eVd8GrmMQDPdQVWdU1YqqWrFs2bKWHy9JaqNNEPwlcEWSs5qjgbXAX7RYbzWwf5L9kuwMHM2gr6Jh5wNPAkiylEFT0fVti5ckbb8Zm4aSBPgK8Fjgtxg09/y3qvo/s224qu5MciLwRWAJcGZVrUtyGrCmqlY2856WZD1wF/BfquoH27VHkqRtMmMQVFUlOb+qHsPWv+ZnVVWrgFVTpp06vH3g5OYlSZqANk1Dlyb5rc4rkSRNRJurhp4EvDTJRuA2Bs1DVVWP6rIwSdJ4tAmCQzuvQpI0MW0eTPOdJI8GHs/gPoCvVtXlnVcmSRqLNncWnwqcDTwQWAp8OMlrui5MkjQebZqGjgEOqqqfwd1dTlwOvKnLwiRJ49HmqqGNwC5D4/cB/rWTaiRJY9fmiOB2YF2SixicI3gq8JUk7wGoqld0WJ8kqWNtguCzzWuLS7opRZI0CW2uGjp7HIVIkiajbe+jkqRFyiCQpJ5rHQRJdu2yEEnSZLS5oex3m26ir23G/32S/9l5ZZKksWhzRPBO4A+BHwBU1VX88kH2kqQFrlXTUFXdMGXSXR3UIkmagDb3EdyQ5HeBah45+QqaZiJJ0sLXJgheCrwb2IvBw+YvBF7WZVHql+WnXDDpEqReaxMEqao/7rwSSdJEtDlH8LUkFyY5PsmenVckSRqrWYOgqvYHXgM8Arg8yeeTHNt5ZZKksWh71dA3qupk4GDghwweVCNJWgTa3FC2e5IXJPkC8DXgewwCQZK0CLQ5WXwVcD5wWlV9veN6JElj1iYIHlpV1XklkqSJmDYIkryrql4JrEyyVRBU1eGdViZJGouZjgjOad7fNo5CJEmTMW0QVNXaZvDAqnr38LwkJwFf6rIwSdJ4tLl89AUjph03z3VIkiZkpnMExwB/BOyXZOXQrN1ouqSWJC18M50j2HLPwFLg7UPTbwGu7rIoSdL4zHSO4DvAd4DfGV85kqRxa3Nn8WOTrE5ya5I7ktyV5CfjKE6S1L02J4vfCxwDfAu4L/Ai4K/bbDzJIUmuS7IhySkzLHdkkkqyos12JUnzp22ncxuAJVV1V1V9GHjSbOskWQKcDhwKHAAck+SAEcvtxuCpZ5dtS+GSpPnRJgh+2jyi8sokb0nyKmDXFusdDGyoquur6g7gXOCIEcu9EXgL8LO2RUuS5k+bIPgTYAlwInAbsA/w3Bbr7QUMP/R+UzPtbkkOAvapqs/PtKEkJyRZk2TN5s2bW3y0JKmtWTuda64eAvh/wBu2YdsZtbm7Zyb3At5Ji5vTquoM4AyAFStW2AGeJM2jmW4o+yZDX9xTVdWjZtn2JgZHD1vsDdw4NL4b8EjgkiQAD2LQwd3hVbVmlm1LkubJTEcEz9jOba8G9k+yH/Bd4GgGdyoDUFU3M7hZDYAklwCvNgQkabxmu6FszqrqziQnAl9kcI7hzKpal+Q0YE1VrZx5C5KkcZj1HEGSW/hlE9HOwL2B26pq99nWrapVwKop006dZtknzrY9SdL8a3OyeLfh8STPwmcWS9Ki0eqGsmFVdT7w5A5qkSRNQJumoecMjd4LWMEMVxNJkhaWNg+vf+bQ8J3ARkbfISxJWoDanCN44TgKkSRNRpumof2AlwPLh5evqsO7K0uSNC5tmobOBz4EfA74RbflSJLGrU0Q/Kyq3tN5JZKkiWgTBO9O8jrgQuD2LROr6vLOqpIkjU2bIPhNBl1RP5lfNg0V3ksgSYtCmyB4NvDQ5uEykqRFps2dxVcBe3ZdiCRpMtocEfwq8M9JVnPPcwRePipJi0CbIHhd51VIkiamzZ3FXxpHIZKkyej0eQSSpB2fzyOQpJ7zeQSS1HM+j0CSes7nEUhSz/k8AknquVnPESQ5O8meQ+P3T3Jmt2VJksalzcniR1XVj7eMVNWPgIO6K0mSNE5tguBeSe6/ZSTJA2h3bkGStAC0+UJ/O/C1JJ9icLXQ84H/0WlVkqSxaXOy+CNJ1jC4dyDAc6pqfeeVSZLGolUTT/PF75e/JC1C23xnsSRpcTEIJKnnDAJJ6jmDQJJ6ziCQpJ7rNAiSHJLkuiQbkpwyYv7JSdYnuTrJPyR5SJf1SJK21lkQJFkCnA4cChwAHJPkgCmLXQGsqKpHAZ8C3tJVPZKk0bo8IjgY2FBV11fVHcC5TOm+uqourqqfNqOXAnt3WI8kaYQug2Av4Iah8U3NtOkcD3xh1IwkJyRZk2TN5s2b57FESVKXQZAR00Y+2SzJsQyefPbWUfOr6oyqWlFVK5YtWzaPJUqSuuxFdBOwz9D43sCNUxdK8hTgz4Hfr6rbO6xnUVl+ygWTLkHSItHlEcFqYP8k+yXZGTgaWDm8QJKDgL8BDq+qmzqsRZI0jc6CoKruBE4EvghcC5xXVeuSnJbk8GaxtwL3Az6Z5MokK6fZnCSpI50+YKaqVgGrpkw7dWj4KV1+viRpdt5ZLEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPbfTpAsYp+WnXDDpEiRph+MRgST1nEEgST1nEEhSz3UaBEkOSXJdkg1JThkx/z5JPtHMvyzJ8i7rkSRtrbMgSLIEOB04FDgAOCbJAVMWOx74UVU9DHgn8Fdd1SNJGq3LI4KDgQ1VdX1V3QGcCxwxZZkjgLOb4U8Bf5AkHdYkSZqiy8tH9wJuGBrfBPz2dMtU1Z1JbgYeCHx/eKEkJwAnNKO3Jrmuk4rnz1Km7MMCtVj2A9yXHdVi2Zex7Ee2r83kIdPN6DIIRv2yrzksQ1WdAZwxH0WNQ5I1VbVi0nVsr8WyH+C+7KgWy74s9P3osmloE7DP0PjewI3TLZNkJ2AP4Icd1iRJmqLLIFgN7J9kvyQ7A0cDK6cssxJ4QTN8JPCPVbXVEYEkqTudNQ01bf4nAl8ElgBnVtW6JKcBa6pqJfAh4JwkGxgcCRzdVT1jtmCasWaxWPYD3Jcd1WLZlwW9H/EHuCT1m3cWS1LPGQSS1HMGwTyarUuNhSLJmUluSnLNpGvZXkn2SXJxkmuTrEty0qRrmoskuyT5RpKrmv14w6Rr2l5JliS5IsnnJ13L9kiyMck3k1yZZM2k65kLzxHMk6ZLjX8BnsrgstjVwDFVtX6ihc1BkicAtwIfqapHTrqe7ZHkwcCDq+ryJLsBa4FnLbT/Ls0d97tW1a1J7g18BTipqi6dcGlzluRkYAWwe1U9Y9L1zFWSjcCKqlqwN8Z5RDB/2nSpsSBU1ZdZJPdzVNX3quryZvgW4FoGd7QvKDVwazN67+a1YH/FJdkbOAz44KRrkUEwn0Z1qbHgvnAWs6Z324OAyyZbydw0TSlXAjcBF1XVgtyPxruA/wr8YtKFzIMCLkyytukOZ8ExCOZPq+4yNBlJ7gd8GnhlVf1k0vXMRVXdVVUHMrhL/+AkC7LZLskzgJuqau2ka5knj6uqRzPoafllTdPqgmIQzJ82XWpoApo29U8DH6uqz0y6nu1VVT8GLgEOmXApc/U44PCmbf1c4MlJPjrZkuauqm5s3m8CPsugmXhBMQjmT5suNTRmzUnWDwHXVtU7Jl3PXCVZlmTPZvi+wFOAf55sVXNTVX9WVXtX1XIG/07+saqOnXBZc5Jk1+YiBJLsCjwNWHBX2xkE86Sq7gS2dKlxLXBeVa2bbFVzk+TjwNeB30iyKcnxk65pOzwO+BMGvzqvbF5Pn3RRc/Bg4OIkVzP40XFRVS3oyy4XiV8FvpLkKuAbwAVV9XcTrmmbefmoJPWcRwSS1HMGgST1nEEgST1nEEhSzxkEktRzBoF2aEkuSdL5Q8GTvKLpofRjU6Yf2OZy0yQrkrynGT4uyXtbfu6/TfKpbaz195oeSK9Mste2rj/Ddt/abPetSV6f5NXN9NOSPGU+PkM7ps4eVSlNWpKdmvs72vhPwKFV9e0p0w9k0EPmqplWrqo1wDZ3QdzclXrkNq72x8DbqurDzfi2rj+dlwDLqur2JK/fMrGqTp2n7WsH5RGBtluS5c2v6Q80vygvbO5+vccv+iRLm24FtvxqPj/J55J8O8mJSU5u+qe/NMkDhj7i2CRfS3JNkoOb9XdtnpuwulnniKHtfjLJ54ALR9R6crOda5K8spn2fuChwMokrxpadmfgNOCo5tf3UUkObmq5onn/jWbZJ47qVz/J85rPuirJl6f5210zVPtnkvxdkm8lecuI5V8EPB84NcnHpqx/WZJHDC17SZLHTPe3mrLdlcCuwGVJjpoy76wkRzbDG5P8VQbPRvhGkoe12U/t4KrKl6/tegHLgTuBA5vx84Bjm+FLGPTVDrAU2NgMHwdsAHYDlgE3Ay9t5r2TQedwW9b/QDP8BOCaZvgvhj5jTwbPgti12e4m4AEj6nwM8M1mufsB64CDmnkbgaUj1jkOeO/Q+O7ATs3wU4BPN8NPBD4/dZ3m8/baUuc0f7trhta7HtgD2AX4DrDPiHXOAo4csf6rgDc0ww8G/mWmv9WI7d46NPx64NUjPm8j8OfN8J8O7fOM++lrx355RKD58u2qurIZXsvgC2o2F1fVLVW1mUEQfK6Z/s0p638c7n5Owu5NnztPA07JoFvmSxh8ce7bLH9RVY16nsLjgc9W1W016Nv/M8Dvtdu9u+0BfLL5Ff5O4BGzLP9V4KwkLwaWtNj+P1TVzVX1M2A98JBtqO084HnN8POBTzbDM/2t5uLjQ++/0wxv635qB+I5As2X24eG7wLu2wzfyS+bIHeZYZ1fDI3/gnv+vzm1H5Ri0O33c6vquuEZSX4buG2aGkd1Fb6t3sggwJ6dwfMNLplp4ap6aVPTYcCVSQ6sqh/MsMrUv2Prf6NV9d0kP0jyKOAoBm3+MM3fajvU1OE57Kd2IB4RqGsbGTTJwNxPah4FkOTxwM1VdTODzv1eniTNvINabOfLwLOS/EoGPUU+G/inWda5hUHz1RZ7AN9tho+b7QOT/FpVXVaDE67f555dlXfhXAYPfNmjqr7ZTJvL32omRw29f73Z5rj3U/PIIFDX3gb8xyRfY3COYC5+1Kz/fmBLT6hvZPC4xqubZpo3zraRGjyy8iwGvUReBnywqq6YZbWLgQO2nCwG3gL8ZZKv0q4J5K0ZPNj8GgZBdFWLdbbHpxh07Xze0LRt/lvN4j5JLgNOYnBeAsa/n5pH9j4qqbUsgge1a2seEUhSz3lEIEk95xGBJPWcQSBJPWcQSFLPGQSS1HMGgST13P8HdxrBW9Jeo34AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "cumprob = np.cumsum(counttails / N)\n", "print('cumprob:', cumprob)\n", "plt.bar(range(0, 6), cumprob, width=1)\n", "plt.xlabel('number of tails in five flips')\n", "plt.ylabel('cumulative probability');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Back to Exercise 2\n", "\n", "Answers to Exercise 3" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Theoretical probabilities:\n", "Exact probability: [0.03125 0.15625 0.3125 0.3125 0.15625 0.03125]\n", "Probability with 1000 trials: [0.031 0.16 0.307 0.334 0.139 0.029]\n", "Probability with 10000 trials: [0.0315 0.1579 0.3083 0.3098 0.1634 0.0291]\n", "Probability with 100000 trials: [0.03093 0.15591 0.31454 0.31266 0.15484 0.03112]\n" ] } ], "source": [ "from scipy.special import comb\n", "print('Theoretical probabilities:')\n", "pexact = np.empty(6)\n", "for k in range(6):\n", " pexact[k] = comb(5, k) * 0.5**k * 0.5**(5 - k)\n", "print('Exact probability:', pexact)\n", "for N in (1000, 10000, 100000):\n", " tails = np.sum(rnd.randint(0, 1+1, (5, N)), axis=0)\n", " counttails = np.zeros(6)\n", " for i in range(6):\n", " counttails[i] = np.count_nonzero(tails==i)\n", " print(f'Probability with {N} trials: {counttails / N}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Back to Exercise 3\n", "\n", "Answers to Exercise 4" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAXOklEQVR4nO3df5AndX3n8eeLRSCKIMomZfjhwrlR16CLjhADkhgQMfEgVZIA0Qg56yijHFHKSq2XHMRNqg71/HEaohBdQRNFxXi1dxCRE1BQ0d0FBBclrCvC3nrHGgg/Tn7cwvv+6J7w9bu9M9+Z3d7vzPh8VE1N//h0f989C/Oa/nT3p1NVSJI0bLdxFyBJmpsMCElSJwNCktTJgJAkdTIgJEmddh93ATvL/vvvX0uWLBl3GZI0r6xbt+4nVbW4a92CCYglS5awdu3acZchSfNKkh9tb51dTJKkTgaEJKmTASFJ6mRASJI6GRCSpE4GhCSpkwEhSepkQEiSOhkQkqROC+ZJakmaiSUrLh93CTvNnef/Ti/79QxCktTJgJAkdTIgJEmdDAhJUicDQpLUyYCQJHUyICRJnQwISVInA0KS1MmAkCR16jUgkpyQ5PYkG5Ks6Fh/TpLbktyS5CtJnjOw7vEkN7dfq/usU5K0rd7GYkqyCLgAeBWwCViTZHVV3TbQ7CZgoqp+muSPgfcAp7TrHq6q5X3VJ0maWp9nEEcAG6pqY1U9BlwKnDTYoKquqaqftrM3AAf2WI8kaQb6HM31AODugflNwJFTtH8T8I8D83slWQtsBc6vqv82vEGSM4EzAQ4++OAdLljS9BbSKKiaWp8BkY5l1dkweQMwAfzGwOKDq2pzkkOBq5PcWlU/+JmdVV0EXAQwMTHRuW9J0uz02cW0CThoYP5AYPNwoyTHAX8GnFhVj04ur6rN7feNwLXA4T3WKkka0mdArAGWJjkkyR7AqcDP3I2U5HDgQppwuGdg+X5J9myn9weOAgYvbkuSetZbF1NVbU1yFnAlsAhYVVXrk6wE1lbVauC9wN7A55MA3FVVJwIvAC5M8gRNiJ0/dPeTJKlnvb5ytKquAK4YWnbuwPRx29nuG8BhfdYmSZqaT1JLkjoZEJKkTgaEJKmTASFJ6mRASJI6GRCSpE4GhCSpkwEhSepkQEiSOhkQkqROBoQkqZMBIUnqZEBIkjoZEJKkTgaEJKmTASFJ6mRASJI6GRCSpE4GhCSpkwEhSepkQEiSOhkQkqROBoQkqZMBIUnqZEBIkjoZEJKkTgaEJKmTASFJ6mRASJI6GRCSpE69BkSSE5LcnmRDkhUd689JcluSW5J8JclzBtadnuSO9uv0PuuUJG2rt4BIsgi4AHgNsAw4LcmyoWY3ARNV9SLgMuA97bbPBM4DjgSOAM5Lsl9ftUqSttXnGcQRwIaq2lhVjwGXAicNNqiqa6rqp+3sDcCB7fSrgauq6t6qug+4Cjihx1olSUP6DIgDgLsH5je1y7bnTcA/zmTbJGcmWZtk7ZYtW3awXEnSoN173Hc6llVnw+QNwATwGzPZtqouAi4CmJiY6Ny3NBcsWXH5uEuQZqzPM4hNwEED8wcCm4cbJTkO+DPgxKp6dCbbSpL602dArAGWJjkkyR7AqcDqwQZJDgcupAmHewZWXQkcn2S/9uL08e0ySdIu0lsXU1VtTXIWzS/2RcCqqlqfZCWwtqpWA+8F9gY+nwTgrqo6saruTfKXNCEDsLKq7u2rVknStvq8BkFVXQFcMbTs3IHp46bYdhWwqr/qJElT8UlqSVInA0KS1MmAkCR1MiAkSZ0MCElSJwNCktTJgJAkdTIgJEmdDAhJUicDQpLUyYCQJHUyICRJnQwISVInA0KS1MmAkCR1GikgkizquxBJ0twy6hnEhiTvTbKs12okSXPGqAHxIuCfgI8luSHJmUn26bEuSdKYjRQQVfVgVf1tVf068KfAecCPk1yS5Lm9VihJGouRr0EkOTHJF4H/CrwPOBT47wy9c1qStDDsPmK7O4BrgPdW1TcGll+W5JidX5YkadxGDYg3VtX1gwuSHFVVX6+qs3uoS5I0ZqNepP5Qx7IP78xCJElzy5RnEEleDvw6sDjJOQOr9gF8NkKSFrDpupj2APZu2z19YPkDwMl9FSVJGr8pA6Kqvgp8NcnFVfWjXVSTJGkOmK6L6YNV9Tbgr5PU8PqqOrG3yiRJYzVdF9On2u//pe9CJElzy3RdTOva71/dNeVIkuaK6bqYbgW26VqaVFUvmmb7E2ievF4EfKyqzh9afwzwQZqxnk6tqssG1j0O3NrO3mV31s+fJSsuH3cJ0s+16bqYXjvbHbdDhF8AvArYBKxJsrqqbhtodhdwBvCOjl08XFXLZ/v5kqQdM10X047cuXQEsKGqNgIkuRQ4CfjXgKiqO9t1T+zA50iSejDlk9RJrm+/P5jkgeHv0+z7AODugflN7bJR7ZVkbTu8+O9up74z2zZrt2zZMoNdS5KmM90ZxNHt96dP1W470rXLGWx/cFVtTnIocHWSW6vqB0P1XQRcBDAxMTGTfUuSpjHqYH0keQlwNM0v+eur6qZpNtkEHDQwfyCwedTPq6rN7feNSa4FDgd+MOVGkqSdZtT3QZwLXAI8C9gfuDjJn0+z2RpgaZJDkuwBnAqsHvHz9kuyZzu9P3AUA9cuJEn9G/UM4jTg8Kp6BCDJ+cCNwF9tb4Oq2prkLOBKmttcV1XV+iQrgbVVtTrJy4AvAvsB/zbJu6rqhcALgAvbi9e7AecP3f0kSerZqAFxJ7AX8Eg7vycjdPdU1RUMvXGuqs4dmF5D0/U0vN03gMNGrE2S1IPpHpT7MM01h0eB9UmuaudfBVw/1baSpPltujOIte33dTRdQZOu7aUaSdKcMd1trpfsqkIkSXPLSNcgkiwF/jOwjOZaBABVdWhPdUmSxmzUd1J/AvgIsBV4JfBJnhwKXJK0AI0aEL9QVV8BUlU/qqq/AH6rv7IkSeM26m2ujyTZDbijfbbhfwG/2F9ZkqRxG/UM4m3AU4GzgZcCfwic3ldRkqTxG+kMon2gjfYs4uyqerDXqiRJYzfqWEwT7dvlbgFuTfKdJC/ttzRJ0jiNeg1iFfCWqroOIMnRNHc2TfnKUUnS/DXqNYgHJ8MBoKquB+xmkqQFbLqxmF7STn47yYXAZ2jGYjoFh9uQpAVtui6m9w3Nnzcw7RvcJGkBm24splfuqkIkSXPLqHcx7Zvk/UnWtl/vS7Jv38VJksZn1IvUq2guSv9++/UAzV1MkqQFatTbXP9NVb1uYP5dSW7uoyBJ0tww6hnEw+2zDwAkOQp4uJ+SJElzwahnEG8GPjlw3eE+HItJkha0aQOiHX/peVX14iT7AFTVA71XJkkaq2m7mKrqCeCsdvoBw0GSfj6Meg3iqiTvSHJQkmdOfvVamSRprEa9BvHvaJ6cfsvQct9JLUkL1KgBsYwmHI6mCYrrgI/2VZQkafxGDYhLaB6O+1A7f1q77Pf7KEqSNH6jBsTzqurFA/PXJPlOHwVJkuaGUS9S35Tk1yZnkhwJfL2fkiRJc8GoZxBHAm9Mclc7fzDwvfY1pFVVvllOkhaYUQPihF6rkCTNOSN1MVXVj6b62t52SU5IcnuSDUlWdKw/JsmNSbYmOXlo3elJ7mi/HNZDknaxUa9BzFiSRcAFwGtobpM9LcmyoWZ3AWcAnx7a9pk0b687EjgCOC/Jfn3VKknaVm8BQfOLfUNVbayqx4BLgZMGG1TVnVV1C/DE0LavBq6qqnur6j7gKuzmkqRdqs+AOAC4e2B+U7tsp22b5MzJt9xt2bJl1oVKkrbVZ0CkY1ntzG2r6qKqmqiqicWLF8+oOEnS1PoMiE3AQQPzBwKbd8G2kqSdoM+AWAMsTXJIkj2AU4HVI257JXB8kv3ai9PHt8skSbtIbwFRVVtp3iNxJfA94HNVtT7JyiQnAiR5WZJNwO8BFyZZ3257L/CXNCGzBljZLpMk7SKjPig3K1V1BXDF0LJzB6bX0HQfdW27CljVZ30L0ZIVl4+7BEkLRJ9dTJKkecyAkCR1MiAkSZ0MCElSJwNCktTJgJAkdTIgJEmdDAhJUicDQpLUyYCQJHUyICRJnQwISVInA0KS1MmAkCR1MiAkSZ0MCElSJwNCktTJgJAkdTIgJEmdDAhJUicDQpLUyYCQJHUyICRJnQwISVInA0KS1MmAkCR1MiAkSZ0MCElSJwNCktSp14BIckKS25NsSLKiY/2eST7brv9WkiXt8iVJHk5yc/v10T7rlCRta/e+dpxkEXAB8CpgE7Amyeqqum2g2ZuA+6rquUlOBd4NnNKu+0FVLe+rPknS1Po8gzgC2FBVG6vqMeBS4KShNicBl7TTlwHHJkmPNUmSRtRnQBwA3D0wv6ld1tmmqrYC9wPPatcdkuSmJF9N8oquD0hyZpK1SdZu2bJl51YvST/n+gyIrjOBGrHNj4GDq+pw4Bzg00n22aZh1UVVNVFVE4sXL97hgiVJT+ozIDYBBw3MHwhs3l6bJLsD+wL3VtWjVfXPAFW1DvgB8Cs91ipJGtJnQKwBliY5JMkewKnA6qE2q4HT2+mTgaurqpIsbi9yk+RQYCmwscdaJUlDeruLqaq2JjkLuBJYBKyqqvVJVgJrq2o18HHgU0k2APfShAjAMcDKJFuBx4E3V9W9fdUqSdpWbwEBUFVXAFcMLTt3YPoR4Pc6tvsC8IU+a5MkTc0nqSVJnQwISVInA0KS1MmAkCR1MiAkSZ16vYtpPlmy4vJxlyBJc4pnEJKkTgaEJKmTASFJ6mRASJI6GRCSpE4GhCSpkwEhSepkQEiSOhkQkqROBoQkqZMBIUnqZEBIkjoZEJKkTgaEJKmTASFJ6mRASJI6GRCSpE4GhCSpkwEhSepkQEiSOhkQkqROBoQkqZMBIUnqZEBIkjr1GhBJTkhye5INSVZ0rN8zyWfb9d9KsmRg3Tvb5bcneXWfdUqSttVbQCRZBFwAvAZYBpyWZNlQszcB91XVc4EPAO9ut10GnAq8EDgB+Jt2f5KkXaTPM4gjgA1VtbGqHgMuBU4aanMScEk7fRlwbJK0yy+tqker6ofAhnZ/kqRdZPce930AcPfA/CbgyO21qaqtSe4HntUuv2Fo2wOGPyDJmcCZ7exDSW7fOaX3Zn/gJ+MuYidZKMeyUI4DPJa5qvdjybt3aPPnbG9FnwGRjmU1YptRtqWqLgIumnlp45FkbVVNjLuOnWGhHMtCOQ7wWOaq+XwsfXYxbQIOGpg/ENi8vTZJdgf2Be4dcVtJUo/6DIg1wNIkhyTZg+ai8+qhNquB09vpk4Grq6ra5ae2dzkdAiwFvt1jrZKkIb11MbXXFM4CrgQWAauqan2SlcDaqloNfBz4VJINNGcOp7bbrk/yOeA2YCvw1qp6vK9ad6F50x02goVyLAvlOMBjmavm7bGk+YNdkqSf5ZPUkqROBoQkqZMBsQskWZXkniTfHXctOyLJQUmuSfK9JOuT/Mm4a5qtJHsl+XaS77TH8q5x17SjkixKclOS/zHuWnZEkjuT3Jrk5iRrx13PbCV5RpLLkny//X/m5eOuaaa8BrELJDkGeAj4ZFX96rjrma0kzwaeXVU3Jnk6sA743aq6bcylzVj7xP7TquqhJE8Brgf+pKpumGbTOSvJOcAEsE9VvXbc9cxWkjuBiaqa1w/KJbkEuK6qPtbeyfnUqvqXcdc1E55B7AJV9TWau7Tmtar6cVXd2E4/CHyPjifc54NqPNTOPqX9mrd/LSU5EPgd4GPjrkWQZB/gGJo7Namqx+ZbOIABoVlqR949HPjWeCuZvbZL5mbgHuCqqpq3xwJ8EPhT4IlxF7ITFPDlJOva4XTmo0OBLcAn2m6/jyV52riLmikDQjOWZG/gC8DbquqBcdczW1X1eFUtp3lS/4gk87L7L8lrgXuqat24a9lJjqqql9CMBP3Wtot2vtkdeAnwkao6HPi/wDavPJjrDAjNSNtf/wXg76vqH8Zdz87QnvpfSzO0/Hx0FHBi23d/KfBbSf5uvCXNXlVtbr/fA3yR+TmS8yZg08BZ6WU0gTGvGBAaWXth9+PA96rq/eOuZ0ckWZzkGe30LwDHAd8fb1WzU1XvrKoDq2oJzWgEV1fVG8Zc1qwkeVp7AwRtl8zxwLy7+6+q/jdwd5LntYuOpRkZYl7pczRXtZJ8BvhNYP8km4Dzqurj461qVo4C/hC4te27B/iPVXXFGGuarWcDl7QvotoN+FxVzevbQxeIXwK+2Pwtwu7Ap6vqS+Mtadb+A/D37R1MG4E/GnM9M+ZtrpKkTnYxSZI6GRCSpE4GhCSpkwEhSepkQEiSOhkQ6lWSJdsbxTbJyiTHTbP9XyR5x06o44wkvzxCu3+tN8lEkg/thM+e9jiH2i9P8ts7+rlD+1yS5A8G5s9I8tc78zO08PgchMamqs7dhR93Bs0DV5tH3aCq1gI7PNz0LI5zOc2orDvz+ZIlwB8An57JRkkWLZDX/WoWPIPQrrAoyd+27134cvvkMkkuTnJyO/3b7bj51yf50NA7DZYluTbJxiRnTy5M8ob2nQ43J7mwHXxvUbvf77bvFHh7+xkTNA8t3Tz5+QP7eWn7XohvAm8dWP6bk3Uk2TvJJ9p93pLkde3y45N8M8mNST7fjlPF0P4Hj/POJO9q29+a5PlDbfcAVgKntLWe0rZ7Rhr/nOSNbdtPJTkuzbstJmu7KckrO/4Nzgde0e7z7e2yX07ypSR3JHnPQA0PtWc93wJenuTYdr+3pnm3yZ5JjkjyD237k5I8nGSPtpaN7fKzk9zW/rwu3d5/HJq7DAjtCkuBC6rqhcC/AK8bXJlkL+BC4DVVdTSweGj75wOvphmT57wkT0nyAuAUmoHdlgOPA6+n+ev7gKr61ao6DPhEVV1Gcybw+qpaXlUPD+3/E8DZVTXVC13+E3B/VR1WVS8Crk6yP/DnwHHt4HJrgXNG+Hn8pG3/EeBnus+q6jHgXOCzba2fBb5O8xT7C2meyH1F2/zXgBtoQ6093tNonhDfa+gzV9C8m2B5VX2gXbac5md4GE0gHdQufxrw3ao6sj2mi4FT2v3vDvwxcCPNaL609XwXeBlwJE+O8LsCOLz9eb15hJ+L5hgDQrvCD6tqcmiOdTTdHYOeD2ysqh+2858ZWn95VT3avkDmHprhGI4FXgqsaYf9OJZmiOWNwKFJPpzkBGDK0WaT7As8o6q+2i761HaaHgdcMDlTVffR/IJeBny9reF04DlTfV5rcpDDrp9Fl+to3i1wDE2oHJbkAODe9p0WR0/WXVXfB34E/MoI+/1KVd1fVY/QjBM0WfvjNAMyAjyP5t/vn9r5S4BjqmorsKEN6iOA97f1vaKtF+AWmrO2NwBbR6hHc4wBoV3h0YHpx9n22ldmsX2AS9q/iJdX1fOq6i/aX9wvphmd9a1M/wKdMNqLgrraheY9EpM1LKuqN42wr8nj6fpZdPkazS/eV9Ac1xbgZJ78RTzdz2+6OoZreWTgusNU+76OZkju/wf8T5qgOrqtF5oXGF1AE+TrknjNc54xIDQXfJ/mr/4l7fwpI2zzFeDkJL8IkOSZSZ7TdvvsVlVfoOkWmhxi+UHg6cM7aYf6vj/J0e2i12/n874MnDU5k2Q/mu6do5I8t1321CSj/OU+nZ+ptaruBvYHllbVRprXo76DJwPia5N1t59/MHD7VPucge8DSyaPkWawxsmzra8BbwO+WVVbgGfRnA2uT7IbcFBVXUPzIqNnANtcn9HcZkBo7NprAm8BvpTkeuD/APdPs81tNP3/X05yC3AVzQitBwDXtl0+FwPvbDe5GPho10VqmlE2L2gvUg9fn5j0V8B+7cXv7wCvbH8pngF8pq3hBppfkDvqGpoL8zcnmQzLbwGT3TzX0Rzn9e3839DcCHAr8FngjKoaPDuAprtna3sx/u2MqO1++iPg8+3+nwA+OlDTL/HkGcMtwC3VjAC6CPi7dpubgA/Mx1du/rxzNFfNCUn2rqqHkoSmW+KOgYupksbAMwjNFf++/at/PbAvzV1NksbIMwhJUifPICRJnQwISVInA0KS1MmAkCR1MiAkSZ3+P1RI1U3GbT8rAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "dice = rnd.randint(1, 6 + 1, (2, 1000))\n", "highest_dice = np.max(dice, 0)\n", "outcome = np.zeros(6)\n", "for i in range(6):\n", " outcome[i] = np.sum(highest_dice == i + 1) / 1000\n", "plt.bar(np.arange(1, 7), height=outcome, width=1)\n", "plt.xlabel('highest dice in two throws')\n", "plt.ylabel('probability');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Back to Exercise 4\n", "\n", "Answers to Exercise 5" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Outcome for 100 throws: [0.01 0.03 0.19 0.13 0.27 0.37]\n", "Outcome for 1000 throws: [0.035 0.076 0.142 0.189 0.262 0.296]\n", "Outcome for 10000 throws: [0.0258 0.0862 0.1323 0.1924 0.248 0.3153]\n", "Exact probabilities: [0.02777778 0.08333333 0.13888889 0.19444444 0.25 0.30555556]\n" ] } ], "source": [ "for N in [100, 1000, 10000]:\n", " dice = rnd.randint(1, 6 + 1, (2, N))\n", " highest_dice = np.max(dice, axis=0)\n", " outcome = np.zeros(6)\n", " for i in range(6):\n", " outcome[i] = np.sum(highest_dice == i + 1) / N\n", " print('Outcome for', N, 'throws: ', outcome)\n", "# Exact values\n", "exact = np.zeros(6)\n", "for i, j in enumerate(range(1, 12, 2)):\n", " exact[i] = j / 36\n", "print('Exact probabilities: ',exact)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Back to Exercise 5\n", "\n", "Answers to Exercise 6" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "polled for A: 522\n", "The Dog will predict the wrong winner\n" ] } ], "source": [ "rnd.seed(2)\n", "poll = rnd.choice([0, 1], 1000, p=[0.49, 0.51]) # A=0, B=1\n", "polled_for_A = np.count_nonzero(poll == 0)\n", "print('polled for A:', polled_for_A)\n", "if polled_for_A > 500: \n", " print('The Dog will predict the wrong winner')\n", "else:\n", " print('The Dog will predict the correct winner')" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1000 polls of 1000 people\n", "Probability that The Dog predicts candidate A to win: 0.262\n" ] } ], "source": [ "Awins = 0\n", "Bwins = 0\n", "for i in range(1000):\n", " poll = rnd.choice([0, 1], 1000, p=[0.49, 0.51])\n", " polled_for_A = np.count_nonzero(poll == 0)\n", " if polled_for_A > 500: \n", " Awins += 1\n", " else:\n", " Bwins += 1\n", "print('1000 polls of 1000 people')\n", "print('Probability that The Dog predicts candidate A to win:', Awins / 1000)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1000 polls of 5000 people\n", "Probability that The Dog predicts candidate A to win: 0.071\n" ] } ], "source": [ "Awins = 0\n", "Bwins = 0\n", "for i in range(1000):\n", " poll = rnd.choice([0, 1], 5000, p=[0.49, 0.51])\n", " polled_for_A = np.count_nonzero(poll == 0)\n", " if polled_for_A > 2500: \n", " Awins += 1\n", " else:\n", " Bwins += 1\n", "print('1000 polls of 5000 people')\n", "print('Probability that The Dog predicts candidate A to win:', Awins / 1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Back to Exercise 6" ] } ], "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" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": false, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false } }, "nbformat": 4, "nbformat_minor": 1 }