notebooks/12_random_variables.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7e5d3951",
   "metadata": {},
   "source": [
    "### Average"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6c12b32f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.5\n",
      "7.0\n",
      "10.0\n"
     ]
    }
   ],
   "source": [
    "def average(numbers):\n",
    "    return sum(numbers) / len(numbers)\n",
    "\n",
    "\n",
    "print(average([1, 2, 3, 4, 5, 6]))\n",
    "print(average([7, 7, 7]))\n",
    "print(average([9.5, 10.5]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8c1ee4af",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.5\n",
      "7\n",
      "10.0\n"
     ]
    }
   ],
   "source": [
    "from statistics import mean\n",
    "\n",
    "\n",
    "print(mean([1, 2, 3, 4, 5, 6]))\n",
    "print(mean([7, 7, 7]))\n",
    "print(mean([9.5, 10.5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82aa36f4",
   "metadata": {},
   "source": [
    "**Problem.** HR management of a company has decided to optimize the average performance of employees by applying the following strategy: while there is a person whose efficiency is less than the average efficiency, fire this person. What will be a result of such strategy?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a7163815",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean of [6, 9, 5, 2, 4, 3] is 4.83\n",
      "Remove 2\n",
      "The mean of [6, 9, 5, 4, 3] is 5.40\n",
      "Remove 3\n",
      "The mean of [6, 9, 5, 4] is 6.00\n",
      "Remove 4\n",
      "The mean of [6, 9, 5] is 6.67\n",
      "Remove 5\n",
      "The mean of [6, 9] is 7.50\n",
      "Remove 6\n"
     ]
    }
   ],
   "source": [
    "from statistics import mean\n",
    "\n",
    "\n",
    "def improve_average(numbers):\n",
    "    minimal, average = min(numbers), mean(numbers)\n",
    "    while minimal < average:\n",
    "        print(f'The mean of {numbers} is {average:.2f}')\n",
    "        print(f'Remove {minimal}')\n",
    "        numbers.remove(minimal)\n",
    "        minimal, average = min(numbers), mean(numbers)\n",
    "\n",
    "\n",
    "improve_average([6, 9, 5, 2, 4, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3dd2265",
   "metadata": {},
   "source": [
    "Average of 1000 throws of a dice:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d0beb050",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.50433\n"
     ]
    }
   ],
   "source": [
    "from random import randint, seed\n",
    "from statistics import mean\n",
    "\n",
    "seed(19)\n",
    "print(mean([randint(1, 6) for _ in range(10 ** 5)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1cc364b",
   "metadata": {},
   "source": [
    "### Birthday Problem\n",
    "\n",
    "Denote by $p(n)$ the probability that \n",
    "among $n$ randomly selected persons there are two with the same birthday.\n",
    "Instead of computing this probability directly, we compute the probability $q(n)$ that among $n$ randomly selected people no two share the same birthday. Clearly, $p(n)=1-q(n)$. One can compute $q(n)$ as follows:\n",
    "$$q(n)=\\frac{365}{365} \\cdot \\frac{364}{365} \\cdot \\frac{363}{365} \\dotsb \\frac{(365-n+1)}{365} \\, .$$\n",
    "Indeed, the first person may have a birthday on any of $365$ days, \n",
    "the second one may have a birthday on all days except for the birthday \n",
    "of the first person and so on. Computing $q(n)$ by hand is not that easy.\n",
    "Instead, let us compute $q(n)$ using Python.\n",
    "\n",
    "The code below computes the value of $p(n)$\n",
    "for all $n=1,2,\\dotsc, 365$. It then plots the function $p(n)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7a575f0b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The probability for 22 people is 0.4756953076625502\n",
      "The probability for 23 people is 0.5072972343239854\n",
      "The probability for 35 people is 0.8143832388747152\n",
      "The probability for 100 people is 0.9999996927510721\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAX80lEQVR4nO3de4xc5X3G8e9vr76s77s2xrf1DbCTUHAWQoRLkwCJoZWdpmkFUhqKSFy1uEqVtClRKpJQVVUSNZEiUVKjAEmahjhJEyzVKaEEhAoBvAZjfMEwNt4btvfqXc+Ovddf/5hjMqx32Vl7ds+cc56PtPLMmcPOw2v78dl33nOOuTsiIhJ9JWEHEBGRwlChi4jEhApdRCQmVOgiIjGhQhcRiYmysN64urraa2trw3p7EZFI2rNnT7u714z2WmiFXltbS319fVhvLyISSWbWMNZrmnIREYkJFbqISEyo0EVEYkKFLiISEyp0EZGYGLfQzewhM2s1s/1jvG5m9h0zS5nZPjPbUPiYIiIynnyO0B8BNr3L67cAa4OvrcADFx9LREQmatx16O7+jJnVvssuW4AfePY6vM+b2VwzW+zuxwsVMgkGh4bp7O2nLd3HqcwAZ/qHODMwxNmBIc4ODtM3MIQ7OM6ww7B79rmf/1xEituN6xbxe8vmFvz7FuLEoiVAU87z5mDbeYVuZlvJHsWzfPnyArx1NJ0+O8DzRzt54WgHh0+eJtWa5kTPWQrVxWaF+T4iMjkWzp5WtIWeN3ffDmwHqKurS9ShZP/gMP9z4AQ7djfx26MdDA07lWUlrF1UxXWrFrBs/gxqqiqorqpk3swKZlSUMq28lOnlpVSWl1BZVkqJQYkZJWaYZYu7xAyDnG1qc5GkKkShtwDLcp4vDbYJ2SmQX+5t4VtPvE5T5xmWzpvO1htWccPaGjasmEtlWWnYEUUkJgpR6DuBbWb2KPABoFvz51mtPWf54s/38fThNtYvns337ngPH758ISUlOooWkcIbt9DN7MfAh4BqM2sGvgKUA7j7d4FdwK1ACsgAd05W2Cg5+FYPdzz8IqfPDvC1ze/hz69boSIXkUmVzyqX28d53YG7C5YoBuqPdXLnw7upmlbGY3dv5PJLZoUdSUQSILTL58bV0bY0n/lBPTWzKvnRZz/A4jnTw44kIgmhU/8LqDszwJ2P7KbUjEfuvFZlLiJTSkfoBfSVnftp6TrDT/7ygyxfMCPsOCKSMDpCL5Bdrx7nl3vfYttH1vD+FfPCjiMiCaRCL4B03yD3PnaA9y2Zw90fXhN2HBFJKE25FMD2Z47Snu7jwU+/n/JS/RspIuFQ+1ykkz1nefCZo/zhlYu5ermmWkQkPCr0i/TA00cYGBrmix+7POwoIpJwKvSL0J0ZYEd9E5uvupQVC2aGHUdEEk6FfhF+9GIDmf4hPrNxVdhRRERU6Beqf3CY7z93jI1rqll/6eyw44iIqNAv1JOHTnKyp4+7Nq4MO4qICKBCv2A/f6mZRbMrueGymrCjiIgAKvQL0p7u4+nDbXz86iWU6pK4IlIkVOgX4LG9bzE47Hxyw9Kwo4iIvE2FfgF+8XIzVy6dw9pFus65iBQPFfoENXdl2N/Swx9duTjsKCIi76BCn6AnDp4E4Ob1l4ScRETknVToE/TrAye5bFEVK6t1ZqiIFBcV+gScyvTz4rFObl6/KOwoIiLnUaFPwG9ea2Vo2PmopltEpAip0Cfg6cNt1Myq5H1L5oQdRUTkPCr0PA0PO8+m2tm4ppoSnUwkIkVIhZ6nwydP09Hbz/VrqsOOIiIyKhV6np5NtQNw/ZoFIScRERmdCj1P/5dqZ3XNTBbPmR52FBGRUanQ89A/OMyLb3ayUdMtIlLEVOh52Nd8ikz/EB9crUIXkeKlQs9DfUMXANfUzgs5iYjI2FToeag/1sWq6pksqKoMO4qIyJhU6ONwd15q7GLDCh2di0hxy6vQzWyTmR02s5SZ3TPK68vN7Ckze9nM9pnZrYWPGo6j7b109vZTp0IXkSI3bqGbWSlwP3ALsB643czWj9jtH4Ed7n41cBvwb4UOGpY9wfx5nebPRaTI5XOEfi2Qcvej7t4PPApsGbGPA7ODx3OAtwoXMVx7jnUxZ3o5q6qrwo4iIvKu8in0JUBTzvPmYFuurwKfMrNmYBfwN6N9IzPbamb1Zlbf1tZ2AXGn3p7GLjYsn6vrt4hI0SvUh6K3A4+4+1LgVuCHZnbe93b37e5e5+51NTU1BXrryZPuG+RIW5qrlmm6RUSKXz6F3gIsy3m+NNiW6y5gB4C7/xaYBkT+LJwDLd24w5VLdblcESl++RT6bmCtma00swqyH3ruHLFPI3AjgJmtI1vo0ZhTeRevtnQD8F5d/1xEImDcQnf3QWAb8DhwiOxqlgNmdp+ZbQ52+wLwWTN7Bfgx8Bfu7pMVeqrsa+5m8Zxp1MzSCUUiUvzK8tnJ3XeR/bAzd9u9OY8PAtcXNlr49rd06+5EIhIZOlN0DD1nBzja3qv5cxGJDBX6GPYH8+fvWzo33CAiInlSoY/h7ULXlIuIRIQKfQwH3+ph8ZxpzJ9ZEXYUEZG8qNDH8NqJ06xbPHv8HUVEioQKfRT9g8OkWtNcccmssKOIiORNhT6KI21pBoedK3SELiIRokIfxWsnegBYpyN0EYkQFfooXjt+morSElZWzww7iohI3lToozh04jRrF1VRVqrhEZHoUGON4tDxHq64RPPnIhItKvQR2tN9tJ3uY91izZ+LSLSo0Ed4/eRpAC5bpEIXkWhRoY9wpDUNwNpFuoeoiESLCn2EVGuaqsoyLpk9LewoIiITokIfIdWWZnXNTMx0U2gRiRYV+gip1jSrF2q6RUSiR4We4/TZAU729LFGhS4iEaRCz3GkrReA1TUqdBGJHhV6jlSwwkVH6CISRSr0HKnWNOWlxor5M8KOIiIyYSr0HKnWNLULZuoaLiISSWquHEfa0ppuEZHIUqEH+gaHaOzMqNBFJLJU6IGGjgxDw64VLiISWSr0gFa4iEjUqdAD5wp9VY3uUiQi0aRCD6Ra0yyZO50ZFWVhRxERuSAq9MCxjl4dnYtIpKnQAw0dGZbrhCIRibC8Ct3MNpnZYTNLmdk9Y+zzZ2Z20MwOmNl/Fjbm5OrODNB9ZoAVC1ToIhJd404Ym1kpcD9wM9AM7Dazne5+MGeftcCXgOvdvcvMFk5W4MnQ2JkBYPl8TbmISHTlc4R+LZBy96Pu3g88CmwZsc9ngfvdvQvA3VsLG3NyNXRmr7KoKRcRibJ8Cn0J0JTzvDnYlusy4DIze9bMnjezTaN9IzPbamb1Zlbf1tZ2YYknwdtH6JpyEZEIK9SHomXAWuBDwO3Ag2Y2d+RO7r7d3evcva6mpqZAb33xGjsyVFdVUFWpJYsiEl35FHoLsCzn+dJgW65mYKe7D7j7m8DrZAs+ErTCRUTiIJ9C3w2sNbOVZlYB3AbsHLHPL8kenWNm1WSnYI4WLubkauxUoYtI9I1b6O4+CGwDHgcOATvc/YCZ3Wdmm4PdHgc6zOwg8BTw9+7eMVmhC6l/cJjj3WdYvkArXEQk2vKaNHb3XcCuEdvuzXnswOeDr0hp7sow7OguRSISeYk/U7RBK1xEJCYSX+hNQaHrCF1Eoi7xhd7QkWF6eSk1syrDjiIiclFU6MGSRTMLO4qIyEVJfKE3dWZYpukWEYmBRBe6u9PYmdFVFkUkFhJd6G2n+zgzMKRCF5FYSHShn1uyqCkXEYmDRBd6Y4eWLIpIfCS60Bs6M5QYLJ2nQheR6Et0oTd29LJ4znQqyhI9DCISE4luMl1lUUTiJPGFrhUuIhIXiS30dN8g7el+XZRLRGIjsYV+7qJcmnIRkbhIbKE3vL1kUTe2EJF4SGyhN3b2AroOuojER2ILvaEjw9wZ5cyZXh52FBGRgkhsoWvJoojEjQpdRCQmElnog0PDtHSd0Rp0EYmVRBb68e6zDA67jtBFJFYSWejnliwu15JFEYmRZBZ6sGRRUy4iEieJLPTGjgwVpSUsmj0t7CgiIgWTzELvzLB0/nRKSyzsKCIiBZPIQm/oyOguRSISO4krdHcPLpurD0RFJF4SV+hdmQHSfYO6MbSIxE7iCr2hI1jhokIXkZhJXKE3BtdB15JFEYmbvArdzDaZ2WEzS5nZPe+y35+YmZtZXeEiFlZjcFKRplxEJG7GLXQzKwXuB24B1gO3m9n6UfabBXwOeKHQIQupoTPDotmVTCsvDTuKiEhB5XOEfi2Qcvej7t4PPApsGWW/fwK+DpwtYL6Ca+zI6C5FIhJL+RT6EqAp53lzsO1tZrYBWObu//1u38jMtppZvZnVt7W1TThsITR09uouRSISSxf9oaiZlQDfAr4w3r7uvt3d69y9rqam5mLfesLODgxxsqdPV1kUkVjKp9BbgGU5z5cG286ZBbwXeNrMjgHXATuL8YPRJq1wEZEYy6fQdwNrzWylmVUAtwE7z73o7t3uXu3ute5eCzwPbHb3+klJfBF+d9lcFbqIxM+4he7ug8A24HHgELDD3Q+Y2X1mtnmyAxbSuTXoKnQRiaOyfHZy913ArhHb7h1j3w9dfKzJ0diZoaqyjPkzK8KOIiJScIk6U7Sho5fl82dgpsvmikj8JKrQs1dZ1HSLiMRTYgp9eNhp6jqj+XMRia3EFPqJnrP0Dw7rpCIRia3EFPq5JYs67V9E4ioxhd6kJYsiEnOJKfSGzl7KSoxL504LO4qIyKRITqF3ZFgybzplpYn5XxaRhElMuzV1ZjTdIiKxlphCb1Chi0jMJaLQu88McCozoJOKRCTWElHojW9fZVFLFkUkvpJR6FqyKCIJkIhCb+jsBdBZoiISa4ko9MaODNVVFVRV5nW1YBGRSEpGoWuFi4gkQCIKvaFDhS4i8Rf7Qu8bHOJ49xlWLNAKFxGJt9gXelPnGYYdaqt1hC4i8Rb7Qj/Wnl3hUqsjdBGJufgXeocKXUSSIRGFPmd6OfNmVoQdRURkUsW+0Bs6MtTqhCIRSYDYF/qb7b1a4SIiiRDrQu8bHOKtU2eorVahi0j8xbrQm7uCJYuachGRBIh1ob+9ZFFH6CKSAPEu9OA66FqyKCJJEO9Cb+9l1rQy5s0oDzuKiMiki3ehd/SysnomZhZ2FBGRSZdXoZvZJjM7bGYpM7tnlNc/b2YHzWyfmT1pZisKH3XijnVoyaKIJMe4hW5mpcD9wC3AeuB2M1s/YreXgTp3vxL4GfCNQgedqP7BYVq6zrBSK1xEJCHyOUK/Fki5+1F37wceBbbk7uDuT7l7Jnj6PLC0sDEnrrkrw7CjI3QRSYx8Cn0J0JTzvDnYNpa7gF+N9oKZbTWzejOrb2tryz/lBXj7oly6bK6IJERBPxQ1s08BdcA3R3vd3be7e52719XU1BTyrc/zZruWLIpIsuRz1+QWYFnO86XBtncws5uALwN/4O59hYl34VKtaebOKGe+rrIoIgmRzxH6bmCtma00swrgNmBn7g5mdjXw78Bmd28tfMyJO9KWZk1NlZYsikhijFvo7j4IbAMeBw4BO9z9gJndZ2abg92+CVQBPzWzvWa2c4xvN2WOtKZZXVMVdgwRkSmTz5QL7r4L2DVi2705j28qcK6L0tXbT0dvP2sWqtBFJDlieabokbY0gApdRBIlloWeas0WuqZcRCRJYlvolWUlLJk3PewoIiJTJp6F3pZmVU0VpSVa4SIiyRHLQj/SlmZ1jU4oEpFkiV2hnx0YornrjD4QFZHEiV2hv3EyjTusXTgr7CgiIlMqdoV+6HgPAOsWq9BFJFniV+gnepheXqrL5opI4sSv0I/3cPkls7TCRUQSJ1aF7u68duK0pltEJJFiVegnes5yKjPAusWzw44iIjLlYlXorx0/DcAVl6jQRSR5YlXoB4MVLldoykVEEihWhX7oeA9L5k5n9rTysKOIiEy5WBX6gbd6WH+ppltEJJliU+idvf282d7LhuXzwo4iIhKK2BT63qYuAK5ePjfcICIiIYlNob/ceIrSEuPKpXPCjiIiEorYFPpLjV1cccksZlTkdZtUEZHYiUWhDw07rzR1a7pFRBItFoWeak2T7hvUB6IikmixKPQ9Dec+EFWhi0hyxaLQnzvSzqLZldQumBF2FBGR0ES+0IeHneeOdHD9mmrMdMlcEUmuyBf6weM9dPb2s3FNddhRRERCFflCf+LgSczg99fWhB1FRCRUkS/0xw+c4JoV86mZVRl2FBGRUEW60FOtaV47cZqPvmdR2FFEREIX6ULfUd9EWYmx+apLw44iIhK6yBZ6pn+Qn+1p5sZ1C1k4a1rYcUREQpdXoZvZJjM7bGYpM7tnlNcrzewnwesvmFltwZOO8PCzx+js7WfrDasn+61ERCJh3EI3s1LgfuAWYD1wu5mtH7HbXUCXu68Bvg18vdBBc+1v6eY7T77BTesW8v4VOjtURATyO0K/Fki5+1F37wceBbaM2GcL8P3g8c+AG22SzvLZUd/EJx54jnkzKviXT1w5GW8hIhJJ+RT6EqAp53lzsG3Ufdx9EOgGFoz8Rma21czqzay+ra3tggKvrJ7JTesWsnPb9VqqKCKSY0ovHu7u24HtAHV1dX4h3+Oa2vlcUzu/oLlEROIgnyP0FmBZzvOlwbZR9zGzMmAO0FGIgCIikp98Cn03sNbMVppZBXAbsHPEPjuBO4LHnwR+4+4XdAQuIiIXZtwpF3cfNLNtwONAKfCQux8ws/uAenffCXwP+KGZpYBOsqUvIiJTKK85dHffBewase3enMdngT8tbDQREZmIyJ4pKiIi76RCFxGJCRW6iEhMqNBFRGLCwlpdaGZtQMMF/ufVQHsB40wW5Sws5SycKGQE5RzNCncf9RZtoRX6xTCzenevCzvHeJSzsJSzcKKQEZRzojTlIiISEyp0EZGYiGqhbw87QJ6Us7CUs3CikBGUc0IiOYcuIiLni+oRuoiIjKBCFxGJicgV+ng3rA6TmR0zs1fNbK+Z1Qfb5pvZE2b2RvDrlN8E1cweMrNWM9ufs23UXJb1nWB895nZhhAzftXMWoLx3Gtmt+a89qUg42Ez+9hUZAzed5mZPWVmB83sgJl9LthebOM5Vs6iGlMzm2ZmL5rZK0HOrwXbVwY3nE8FN6CvCLZP+Q3p3yXjI2b2Zs5YXhVsD+X3HAB3j8wX2cv3HgFWARXAK8D6sHPl5DsGVI/Y9g3gnuDxPcDXQ8h1A7AB2D9eLuBW4FeAAdcBL4SY8avA342y7/rg974SWBn8mSidopyLgQ3B41nA60GeYhvPsXIW1ZgG41IVPC4HXgjGaQdwW7D9u8BfBY//Gvhu8Pg24CchZnwE+OQo+4fye+7ukTtCz+eG1cUm9wba3wc+PtUB3P0ZstepzzVWri3ADzzreWCumS0OKeNYtgCPunufu78JpMj+2Zh07n7c3V8KHp8GDpG9p26xjedYOccSypgG45IOnpYHXw58hOwN5+H88ZySG9LnkXEsofyeQ/SmXPK5YXWYHPi1me0xs63BtkXufjx4fAJYFE6084yVq9jGeFvwY+tDOdNVRZEx+HH/arJHbEU7niNyQpGNqZmVmtleoBV4guxPB6c8e8P5kVnyuiH9ZGd093Nj+c/BWH7bzM7dtT60sYxaoRe7je6+AbgFuNvMbsh90bM/jxXdOtFizQU8AKwGrgKOA/8aapocZlYF/Bz4W3fvyX2tmMZzlJxFN6buPuTuV5G9X/G1wBXhJjrfyIxm9l7gS2SzXgPMB/4hvIRZUSv0fG5YHRp3bwl+bQV+QfYP58lzP24Fv7aGl/AdxspVNGPs7ieDv0jDwIP8bgog1IxmVk62JH/k7v8VbC668RwtZ7GOaZDtFPAU8EGy0xTn7qiWmyXUG9LnZNwUTGu5u/cBD1MEYxm1Qs/nhtWhMLOZZjbr3GPgo8B+3nkD7TuAx8JJeJ6xcu0EPh18Un8d0J0zlTClRsw7/jHZ8YRsxtuCFQ8rgbXAi1OUycjeQ/eQu38r56WiGs+xchbbmJpZjZnNDR5PB24mO9//FNkbzsP54zmlN6QfI+NrOf+AG9k5/tyxDOfv0FR9+lqoL7KfIL9Odp7ty2Hnycm1iuwqgVeAA+eykZ3fexJ4A/hfYH4I2X5M9sfrAbLzeXeNlYvsJ/P3B+P7KlAXYsYfBhn2kf1Lsjhn/y8HGQ8Dt0zhWG4kO52yD9gbfN1ahOM5Vs6iGlPgSuDlIM9+4N5g+yqy/6CkgJ8ClcH2acHzVPD6qhAz/iYYy/3Af/C7lTCh/J67u079FxGJi6hNuYiIyBhU6CIiMaFCFxGJCRW6iEhMqNBFRGJChS4iEhMqdBGRmPh/TbFT+9SXGB8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from itertools import accumulate\n",
    "import matplotlib.pyplot as plt\n",
    "from operator import mul\n",
    "\n",
    "num_days = 365\n",
    "\n",
    "factors = [1 - i / num_days for i in range(num_days)]\n",
    "prob_no_collision = list(accumulate(factors, mul))\n",
    "prob_collision = [1 - p for p in prob_no_collision]\n",
    "\n",
    "plt.plot(prob_collision)\n",
    "plt.savefig('birthdays.png')\n",
    "\n",
    "for n in (22, 23, 35, 100):\n",
    "    print(f'The probability for {n} people '\n",
    "          f'is {prob_collision[n - 1]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0bd6261",
   "metadata": {},
   "source": [
    "### Expectation Is Not All\n",
    "\n",
    "Alice and Bob play the following game. Each of them has fair, but unconventional dice. The numbers on the sides of Alice's dice are $2, 2, 2, 2, 3, 3$, whereas the numbers on the sides of Bob's dice are $1, 1, 1, 1, 6, 6$.\n",
    "They throw their dices. The one with the larger number on the dice wins. If Alice and Bob play the game many times, who is likely to win more often?\n",
    "\n",
    "Let us compute the expectations of the corresponding two random variables.\n",
    "Alice gets $2$ with probability $2/3$ and gets $3$ with probability $1/3$. Her expectation is \n",
    "$$\n",
    "2 \\cdot \\frac 23 + 3 \\cdot \\frac 13 = \\frac 73.\n",
    "$$\n",
    "Bob gets $1$ with probability $2/3$ and gets $6$ with probability $1/3$. His expected value is \n",
    "$$\n",
    "1 \\cdot \\frac 23 + 6 \\cdot \\frac 13 = \\frac 83.\n",
    "$$\n",
    "Therefore, Bob has better expected value. Does it mean that he wins more often?\n",
    "\n",
    "Let us compute the probabilities for each player to win directly. \n",
    "The result of the game depends only on the outcome of Bob: indeed, \n",
    "if he throws $1$, he loses, if he throws $6$, he wins. \n",
    "Bob throws $1$ with probability $2/3$. Thus, Bob loses with probability $2/3$. Bob loses more often despite having larger expected value.\n",
    "\n",
    "The following code snippet simulates the game $100\\,000$ times\n",
    "and shows that, indeed, Alice wins in about $2/3$ of all games."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ab2ca5af",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Out of 100000 throws, dice1 won 66868 times\n"
     ]
    }
   ],
   "source": [
    "from random import randint, seed\n",
    "\n",
    "seed(27)\n",
    "\n",
    "dice1 = [2, 2, 2, 2, 3, 3]\n",
    "dice2 = [1, 1, 1, 1, 6, 6]\n",
    "\n",
    "num_rounds = 10 ** 5\n",
    "num_dice1_wins = 0\n",
    "\n",
    "for _ in range(num_rounds):\n",
    "    dice1_result = dice1[randint(0, 5)]\n",
    "    dice2_result = dice2[randint(0, 5)]\n",
    "\n",
    "    if dice1_result > dice2_result:\n",
    "        num_dice1_wins += 1\n",
    "\n",
    "print(f'Out of {num_rounds} throws, dice1 won {num_dice1_wins} times')"
   ]
  }
 ],
 "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}