notes3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 2. Discrete random variables\n",
    "\n",
    "\n",
    " #### [Back to main page](https://petrosyan.page/fall2020math3215)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Example 1\n",
    "\n",
    "Consider the  double coin flip experiment $S=\\{HH, HT, TH, TT\\}$, and define a random variable which counts the number of heads in each outcome. Then \n",
    "$$X(HH)=2,\\;X(HT)=1,\\;X(TH)=1,\\; X(TT)=0.$$\n",
    " \n",
    "Observe that Range$(X)=\\{0,1,3\\}$ and the pmf of the random variable is \n",
    "        \n",
    "$$f(0)=0.25,\\; f(1)=0.5,\\; f(2)=0.15.$$\n",
    "\n",
    "Below we draw the line graph and the histogram of the pmf of this random variable. Red circles represent the range of $X$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAE/CAYAAAC0Fl50AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAASUklEQVR4nO3df7Dld13f8dc7bMBAfhEQhZBAaiMKylSSxmCr0FokRhiYVtsg0qaDpSo6LaPtaGs1QKOD7eCggWqYYgumAYYyDLYJYluJaJNMCGA1QJo0CZOYCCaYbBYCZJN3/zgn8frOze7d7M09u/c+HjN39vz4nnM+97tn3vuc7/2eu9XdAQAA/sIRq14AAAAcakQyAAAMIhkAAAaRDAAAg0gGAIBBJAMAwCCSOaRU1XdW1bWrXsfBqqoXVtUtq14HwFapqmuq6oWrXgdsFpHMSlTVTVX1d+bt3f3R7n7WKtYEwMNbb25X1blV9ftJ0t3P6e6P7Oc5nllVXVW7HsWlwqYQybAOAxzg8GN2s5lEMoeUeZrC8sjFT1XV/6mqu6rqPVX1NWvuf0lVfbKq7qyq/11Vz93Hc39PVV27fJ63VdVlVfXDy/vOrao/qKpfrqovJDmvqr6hqv5XVd1RVbdX1UVVdfxY289U1aeq6s+r6jfWrm25zU9W1eer6raq+seburMADiFrjzRX1RlV9bGq2l1Vn6uqNy83+73ln3dW1Z6qen5VHVFVP1tVn13Oy3dW1XFrnvcfLu+7o6r+zXid86rqfVX1m1W1O8m5y9e+fPnvwm1VdUFVPXbN83VV/VhVXVdVd1fVG5fz/vLlet+7dnt2LpHM4eDvJzkrySlJnpvk3CSpqucleUeSf5rkSUl+PckHq+px8wmq6slJ3pfkZ5bbXpvkO8Zm357khiRPSXJ+kkryi0meluSbk5yU5LzxmFcmeXGSb0jyjUl+ds19X5/kuCQnJnl1krdW1RMP7FsHOCy9JclbuvvYLObje5e3f9fyz+O7++juvjyLmX5ukr+V5K8kOTrJBUlSVc9O8rYsZu1T8xczda2XZTHfj09yUZL7krwuyZOTPD/Jdyf5sfGYs5KcluTMJP8yyYXL1zgpybckecVBfO9sEyKZw8GvdPet3f2FJL+V5K8tb/8nSX69u6/s7vu6+z8n+UoWQ286O8k13f3+7t6b5FeS/OnY5tbu/tXu3tvd93T39d39O939le7+syRvTvKC8ZgLuvvm5drOz18erPcmeUN339vdlyTZk8T51sDh7APLI7R3VtWdWQTseu5N8ler6sndvae7r9jHc74yyZu7+4bu3pPFwYxzlqdOfH+S3+ru3+/uryb5uSQ9Hn95d3+gu+9fzu6ru/uK5Sy/KYsDKHN2v6m7d3f3NUn+OMmHl69/V5JLk3zbxncJ25VI5nCwNma/lMVRhiR5RpKfHAP7pCyO/E5PS3LzA1e6u5PM3z5x89orVfWUqnp3Vf3J8sd4v5nFkYmHe8xnx2vfsQzy9dYOcDh6eXcf/8BXHnqE9gGvzuKna5+pqquq6iX7eM6nZTE/H/DZJLuSfF0eOru/lOSO8fg5u7+xqv5bVf3pcnb/Qh46uz+35vI961w3qxHJHNZuTnL+2oHd3Y/v7ovX2fa2JE9/4EpV1drrS/PoxC8ub3vu8keGP5TFKRhrnbTm8slJbn0E3wfAttLd13X3K7I4fe1NSd5XVU/IQ+dsspibz1hz/eQke7MI1zm7j8rilLm/9HLj+n9I8pkkpy5n97/KQ2c37JdIZpWOrKqvWfN1oJ9KfnuSH6mqb6+FJ1TV91XVMets+9+TfGtVvXz5Oq/N4pzhfTkmi1Mk7qyqE5P8i3W2eW1VPb2qTshiEL/nAL8HgG2nqn6oqr62u+9Pcufy5vuS/FmS+7M49/gBFyd5XVWdUlVHZ3Hk9z3Ln8S9L8lLq+o7lh+me332H7zHJNmdZE9VfVOSH920b4wdRSSzSpdk8WOtB77OO5AHd/fHsjgv+YIkf57k+iw/1LfOtrcn+YEkv5TFj+qeneRjWZzD/HBen+R5Se7KIrLfv842/yXJh7P4wN8NSf7tgXwPANvUWUmuqao9WXyI75zu/vLydInzk/zB8jS5M7P4APa7svjNFzcm+XKSn0iS5TnDP5Hk3VkcVb47yeez79n9U0l+cLnt2+PgBY9QLU7NhJ2lqo7I4pzkV3b37z7C57gpyQ939//YzLUBsL7lkeY7sziV4sZVr4ftzZFkdoyqenFVHb/8FXEPnKO2r09cA7BiVfXSqnr88pzmf5/kj5LctNpVsROIZHaS5yf5f0luT/LSLD6lfc9qlwTAfrwsiw/33Zrk1CxO3fBjcB51TrcAAIDBkWQAABhEMgAADAf6e2m3xFlnndUf+tCHVr0MgEdix/2nBWY2cBh72Jl9SB5Jvv3221e9BAA2yMwGtqNDMpIBAGCVRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGDYUCRX1VlVdW1VXV9VP73O/S+sqruq6pPLr5/b6GMB2FxmNsDB27W/DarqMUnemuRFSW5JclVVfbC7PzU2/Wh3v+QRPhaATWBmA2yOjRxJPiPJ9d19Q3d/Ncm7k7xsg89/MI8F4MCZ2QCbYCORfGKSm9dcv2V52/T8qvrDqrq0qp5zgI8FYHOY2QCbYL+nWySpdW7rcf3jSZ7R3Xuq6uwkH0hy6gYfu3iRqtckeU2SnHzyyRtYFhy811+9+8HLP3/asStcCWwaM5sN2Xt/Z9cR6/2V82ixzw8vG4nkW5KctOb605PcunaD7t695vIlVfW2qnryRh675nEXJrkwSU4//fR1hzJstjd+Ys+Dl0Uy24SZzYbsOqJy9qV35Ma79656KTvCKcfsyiXf+6RVL4MDsJFIvirJqVV1SpI/SXJOkh9cu0FVfX2Sz3V3V9UZWZzGcUeSO/f3WAA2lZnNht14995ct/u+VS8DDkn7jeTu3ltVP57kt5M8Jsk7uvuaqvqR5f2/luT7k/xoVe1Nck+Sc7q7k6z72EfpewHY8cxsgM2xkSPJ6e5Lklwybvu1NZcvSHLBRh8LwKPHzAY4eP7HPQAAGEQyAAAMIhkAAAaRDAAAg0gGAIBBJAMAwCCSAQBgEMkAADCIZAAAGEQyAAAMIhkAAAaRDAAAg0gGAIBBJAMAwCCSAQBgEMkAADCIZAAAGEQyAAAMIhkAAAaRDAAAg0gGAIBBJAMAwCCSAQBgEMkAADCIZAAAGEQyAAAMIhkAAAaRDAAAg0gGAIBBJAMAwCCSAQBgEMkAADCIZAAAGEQyAAAMIhkAAAaRDAAAg0gGAIBBJAMAwCCSAQBgEMkAADCIZAAAGEQyAAAMIhkAAAaRDAAAg0gGAIBBJAMAwCCSAQBgEMkAADCIZAAAGEQyAAAMIhkAAAaRDAAAg0gGAIBBJAMAwCCSAQBgEMkAADCIZAAAGEQyAAAMIhkAAAaRDAAAg0gGAIBBJAMAwCCSAQBgEMkAADCIZAAAGEQyAAAMIhkAAAaRDAAAg0gGAIBBJAMAwCCSAQBgEMkAADCIZAAAGEQyAAAMIhkAAAaRDAAAg0gGAIBBJAMAwCCSAQBgEMkAADCIZAAAGEQyAAAMIhkAAAaRDAAAg0gGAIBBJAMAwCCSAQBgEMkAADCIZAAAGEQyAAAMIhkAAAaRDAAAg0gGAIBBJAMAwCCSAQBgEMkAADCIZAAAGEQyAAAMIhkAAAaRDAAAg0gGAIBBJAMAwCCSAQBgEMkAADCIZAAAGDYUyVX1r6vqq8uvS9e5/61Vdc/ya3dV/cCa+/ZW1ZeX931xMxcPwEOZ2QAHb7+RXFVHJjkvyYuSPDHJC6vqpWOzP0zyzd19VJJ/l+Q3xv3P7e6juvsJB79kAB6OmQ2wOTZyJPncJHd192Xd/cUklyV57doNuvvC7r5pefVdSY7azEUCsGHnxswGOGgbieRnJbl9zfWbkjxtH9v/apLr1lzvJJ+sqi9W1bsOeIUAHAgzG2AT7NrANrXObb3uhlWvS/I9SZ6z5uYzuvsTVfXsJFdX1ZXdfcE6j31Xkr+bJCeccMIGlrU9vf7q3Q9e/vnTjl3hSuDR4T3+qDtsZ/be+zu7jlhv+XD4O+Fx5T2+xQ52f28kkj+T5B+tuf7MJLfNjarq7yX5pSTf193XP3B7d39i+eenqurKJC9O8pCB292vSvKqJDn99NPXHeg7wRs/sefBywKC7ch7/FF32M7sXUdUzr70jtx4997NeDr24W9+3WPz9hc8cdXL2FGOPfII7/EtdMoxu3LJ9z7poJ5jI5H8ziRvq6rvTPLxJC9I8g/WblBVZya5OMlru/vDa27/2iS7uvu25eXnJXnTQa0YgH05rGf2jXfvzXW779vKl9yRnnm0fbwq3uOHj/2ek9zdX0nyxiT/M8mdST7a3R+sqouq6qLlZv8pi+B+y/i1Qc9OckNV3ZPk5iRXdPf5m/1NALBgZgNsjo0cSU53vyHJG8Ztr1xz+Zse5nGXxaemAbaUmQ1w8PyPewAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRDIAAAwiGQAABpEMAACDSAYAgEEkAwDAIJIBAGAQyQAAMIhkAAAYRPKh5OKL932dzWV/bz37HIDDRHX3qtfwEMccc0yfdtppq17G1rr33uSKK3LZsU9JvvKlJMlx9+xOjj02qVrx4rah7mT37tx11LEP3mR/P8rW7vPHPT55yil5wfVXJmeemRx55KpXt2k+8pGP7Lg30GbO7D/6wr25575D79+l7ea4I4/Is47fZX9vIft8ax31mMq3nrD/f1v2NbMPyUiuqg9191mrXseWqvrrSX4nyXEXJnnN4ta7krwo3VetbmHblP299ezzbcvM9n7eEvb51rK/D81I3umq6mPdffqq17FT2N9bzz5nO/F+3nr2+dbaqfvbOckAADCIZAAAGETyoenCVS9gh7G/t559znbi/bz17POttSP3t3OSAQBgcCQZAAAGkXyIqaqzquraqrq+qn561evZzqrqHVX1+ar641WvZSeoqpOq6ner6tNVdU1V/bNVrwkOlpm9tcztrWNmO93ikFJVj0nyf5O8KMktSa5K8oru/tRKF7ZNVdV3JdmT5J3d/S2rXs92V1VPTfLU7v54VR2T5OokL/f+5nBlZm89c3vrmNmOJB9qzkhyfXff0N1fTfLuJC9b8Zq2re7+vSRfWPU6doruvq27P768fHeSTyc5cbWrgoNiZm8xc3vrmNki+VBzYpKb11y/JTvsDcnOUFXPTPJtSa5c7UrgoJjZ7Ag7dWaL5EPLev9/uPNh2Faq6ugk/zXJP+/u3ateDxwEM5ttbyfPbJF8aLklyUlrrj89ya0rWgtsuqo6Mothe1F3v3/V64GDZGazre30mS2SDy1XJTm1qk6pqscmOSfJB1e8JtgUVVVJ/mOST3f3m1e9HtgEZjbblpktkg8p3b03yY8n+e0sTpB/b3dfs9pVbV9VdXGSy5M8q6puqapXr3pN29zfSPKqJH+7qj65/Dp71YuCR8rM3nrm9pba8TPbr4ADAIDBkWQAABhEMgAADCIZAAAGkQwAAINIBgCAQSQDAMAgkgEAYBDJAAAw/H98yaf6lbuQXQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# nbi:hide_in\n",
    "# library\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "plt.rcParams[\"figure.figsize\"] = (12, 5)\n",
    "\n",
    "\n",
    "range_x = np.array([0, 1, 2])\n",
    "pmf_values = np.array([1/4, 1/2, 1/4])\n",
    "\n",
    "fig, [ax1, ax2] = plt.subplots(1,2, num=1, clear=True)\n",
    "\n",
    "ax1.set_ylim(0,0.6) \n",
    "ax1.set_xlim(-0.7, 2.6)\n",
    "ax1.axhline(y=0, color='k')\n",
    "ax1.set_xticks(range_x)\n",
    "ax1.set_yticks(pmf_values)\n",
    "ax1.spines[\"top\"].set_visible(False)  \n",
    "ax1.spines[\"right\"].set_visible(False)\n",
    "ax1.spines[\"bottom\"].set_visible(False)\n",
    "\n",
    "ax2.set_ylim(0, 0.6) \n",
    "ax2.set_xlim(-0.7, 2.6)\n",
    "ax2.axhline(y=0, color='k')\n",
    "ax2.set_xticks(range_x)\n",
    "ax2.set_yticks(pmf_values)\n",
    "ax2.spines[\"top\"].set_visible(False)  \n",
    "ax2.spines[\"right\"].set_visible(False)\n",
    "ax2.spines[\"bottom\"].set_visible(False)\n",
    "\n",
    "# Plotting line graphs using plt.stem with stems removed\n",
    "ax1.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20)\n",
    "markers,stems,base = ax1.stem(range_x, pmf_values, markerfmt=' ', \n",
    "                              linefmt='#039be5', basefmt=\"black\", use_line_collection=True)\n",
    "stems.set_linewidth(3)\n",
    "ax1.set_title(\"Line graph\")\n",
    "\n",
    "\n",
    "# PLotting with plt.bar instead of plt.hist works better when f(x) are knowwn\n",
    "ax2.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20)\n",
    "ax2.bar(range_x, pmf_values, width=1, color='#039be5', edgecolor=\"w\", linewidth=1.3, label=\"Histogran\")\n",
    "ax2.set_title(\"Histogram\")\n",
    "\n",
    "\n",
    "plt.show();\n",
    "                       "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Example 2\n",
    "\n",
    "If Range$(X)=\\{x_1,\\dots,x_k\\}$ and $P(x_1)=\\cdots=P(x_k)=\\frac{1}{k}$ then we say that $X$ has uniform distribution.\n",
    "\n",
    "in the example below, the range of $X$ is \n",
    "\n",
    "$$ \\text{Range}(X)=\\{-3, -6, 5, 9, 2, 1, 3, 8, 7, 10\\}$$\n",
    "\n",
    "with $f(x)=0.1$ for every $x\\in \\text{Range}(X).$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAE/CAYAAABM9qWDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAd7UlEQVR4nO3df5Tdd13n8eeLhBQKFAIUJD+mRg1o1i1Ih4BrRRSqKatG/LGmIghLzKmH+GP3+COyri7LsoKre9Slboyl/ugq2e6haNRK+eHyaxVJwBaaYtlsgGZMpC1EKohJp33vH3PLuUwnc78zc2fu3Hyej3Pume+Pz/d+3/d7733nle/9ztxUFZIkSVJrHjbqAiRJkqRRMAhLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBBWkzLjt5OcTvKBUdcjSYIkX56kkqztzT85yXuS/EOSXxl1fTr/rB11AdKIXA5cAWyqqs+PuhhJ0pz2APcAF5VffKBl4BlhteoS4BOGYEla1S4BbjcEa7kYhDVWknwiyc8mub13WcNvJ3lEkuclmUry00nuSnIqyXcleWGSjyX5TJJX9e7jFcC1wNcn+VySV4/2UUnS+S/J5iQ3Jrk7yaeTvCHJmiS/nOSeJMeBf9k3/neAHwJ+uterXzCq2nX+8tIIjaMXA98GfB74Y+DngHcAXwY8AtgIvAz4LeDtwGXABPDBJAer6o1J7gd2V9XlK1++JLUlyRrgT4A/B14C3A9MAj8MfDvwdcz09Dc/uE1VvSwJwFRV/dxK16w2eEZY4+gNVXWiqj4DvBa4qrf8PuC1VXUfcBB4IvBrVfUPVXUUOApcOpKKJalt24ENwE9V1eer6p+q6n3AvwJ+ta+n/+JIq1RzDMIaRyf6pj/JTHMF+HRV3d+b/kLv56f6xn4BePQy1yZJeqjNwCeranrW8g08tKdLK8YgrHG0uW96Ajg5qkIkSZ2cACYe/LNofU7x0J4urRiDsMbRK5NsSvJ44FXA/xx1QZKkeX2AmdD7uiSP6v2S8zcANwA/1uvp64F9I61SzTEIaxz9AfA24Hjv9p9GW44kaT69y9a+A/gq4E5gCvh+Zn6p+WbgVuBDwI2jqlFtin+aT+MkySeY+WsP7xh1LZIkabx5RliSJElNMghLkiSpSV4aIUmSpCZ5RliSJElNMghLkiSpSbP/sPWK2bFjR731rW8d1e4laaky6gJWkj1b0pibs2eP7IzwPffcM6pdS5IWyJ4t6XzkpRGSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJA4Nwko8leSDJP51jfZLckuRski8k+YHhlylJ6sq+LUnddDkj/OvAD86z/t8DG4ELgL3AgSHUJUlaPPu2JHUwMAhX1RuAO+cZsgs4WDPeCKxL8vRhFShJWhj7tiR1M4xrhJ8A3N43/zng0iHcryRpedi3JQlYO4T7yBzLHphzYHI98N0Aj3/844ewa0nSInTq2/ZsScN08fWnOH2mRrLv6d0b5lw+jCB8D7Ctb/7RwG1zDayqlwAvAZicnBzNkZAkderb9mxJw3T6TJ0zkI7KMC6NuAHY1fst5FcAZ6vq1iHcryRpedi3JYkOZ4STfBLYBDwsyTRwPbAOoKpeDLwa+B7gLHA/8MPLVq0kaSD7tiR1MzAIV9UlA9YX8M+HVpEkaUns25LUjd8sJ0mSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDWpUxBOsiPJHUmOJdk3x/r1Sd6S5MNJPpDka4dfqiSpC3u2JHUzMAgnWQNcA1wJbAOuSrJt1rBXAbdU1aXAS4FfG3ahkqTB7NmS1F2XM8LbgWNVdbyqzgIHgZ2zxmwD3glQVX8DfHmSJw+1UklSF/ZsSepobYcxG4ETffNTwLNnjbkV+G7gfUm2A5cAm4BP9Q9KsgfYAzAxMbHIkiVJ87BnS1o2F19/itNnalHbrr8gQ65m6boE4bmqnn0EXgf8WpJbgI8Afw1MP2SjqgPAAYDJycnFHUVJ0nzs2ZKWzekzxfTuDaMuY2i6BOEpYHPf/CbgZP+AqroXeDlAkgAf790kSSvLni1JHXW5RvgwsDXJliTrgF3Aof4BSR7XWwewG3hPr9FKklaWPVuSOhp4RriqppPsBW4G1gDXVdXRJFf31u8Hvgb4vST3A7cDr1jGmiVJ52DPlqTuulwaQVXdBNw0a9n+vum/BLYOtzRJ0mLYsyWpG79ZTpIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpM6BeEkO5LckeRYkn1zrH9skj9OcmuSo0lePvxSJUld2LMlqZuBQTjJGuAa4EpgG3BVkm2zhr0SuL2qng48D/iVJOuGXKskaQB7tiR11+WM8HbgWFUdr6qzwEFg56wxBTwmSYBHA58BpodaqSSpC3u2JHXUJQhvBE70zU/1lvV7A/A1wEngI8CPV9UDQ6lQkrQQ9mxJ6mhthzGZY1nNmv824BbgW4CvBN6e5L1Vde+X3FGyB9gDMDExsfBqJUmD2LMlzevi609x+szsttDN+gvmajHjq0sQngI2981vYuYsQr+XA6+rqgKOJfk48NXAB/oHVdUB4ADA5OTk4p4BSdJ87NmS5nX6TDG9e8Ooy1gVulwacRjYmmRL75cpdgGHZo25E3g+QJInA08Djg+zUElSJ/ZsSepo4BnhqppOshe4GVgDXFdVR5Nc3Vu/H3gN8DtJPsLMx3I/U1X3LGPdkqQ52LMlqbsul0ZQVTcBN81atr9v+iTwrcMtTZK0GPZsSerGb5aTJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUJIOwJEmSmtQpCCfZkeSOJMeS7Jtj/U8luaV3uy3J/UkeP/xyJUmD2LMlqZuBQTjJGuAa4EpgG3BVkm39Y6rqv1TVM6rqGcDPAu+uqs8sR8GSpHOzZ0tSd13OCG8HjlXV8ao6CxwEds4z/irgTcMoTpK0YPZsSeqoSxDeCJzom5/qLXuIJBcCO4A3L700SdIi2LMlqaO1HcZkjmV1jrHfAfyfc33ElmQPsAdgYmKiU4GSpAWxZ0sNuPj6U5w+c6639vzWXzBXm2hTlyA8BWzum98EnDzH2F3M8xFbVR0ADgBMTk4u7tmTJM3Hni014PSZYnr3hlGXMfa6XBpxGNiaZEuSdcw0zkOzByV5LPBNwB8Nt0RJ0gLYsyWpo4FnhKtqOsle4GZgDXBdVR1NcnVv/f7e0BcBb6uqzy9btZKkedmzJam7VI3m067Jyck6cuTISPYtSUPQ1EV29mxpdVl77UkvjViYOXu23ywnSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNalTEE6yI8kdSY4l2XeOMc9LckuSo0nePdwyJUld2bMlqZu1gwYkWQNcA1wBTAGHkxyqqtv7xjwO+A1gR1XdmeRJy1WwJOnc7NmS1F2XM8LbgWNVdbyqzgIHgZ2zxvwAcGNV3QlQVXcNt0xJUkf2bEnqqEsQ3gic6Juf6i3r91RgfZJ3JflgkpcOq0BJ0oLYsyWpo4GXRgCZY1nNcT+XAc8HHgn8ZZL3V9XHvuSOkj3AHoCJiYmFVytJGsSeLa2Qi68/xekzs99eK2P9BXO91bVQXYLwFLC5b34TcHKOMfdU1eeBzyd5D/B04EuaalUdAA4ATE5OjuaVI0nnN3u2tEJOnymmd28YdRlagi6XRhwGtibZkmQdsAs4NGvMHwHfmGRtkguBZwMfHW6pkqQO7NmS1NHAM8JVNZ1kL3AzsAa4rqqOJrm6t35/VX00yVuBDwMPANdW1W3LWbgk6aHs2ZLUXapG82nX5ORkHTlyZCT7lqQhaOoCPXu29FBrrz3ppRHjY86e7TfLSZIkqUkGYUmSJDXJICxJkqQmGYQlSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJalKnIJxkR5I7khxLsm+O9c9L8tkkt/RuPz/8UiVJXdizJambtYMGJFkDXANcAUwBh5McqqrbZw19b1V9+zLUKEnqyJ4tSd11OSO8HThWVcer6ixwENi5vGVJkhbJni1JHXUJwhuBE33zU71ls319kluT/FmSfzaU6iRJC2XPlqSOBl4aAWSOZTVr/kPAJVX1uSQvBP4Q2PqQO0r2AHsAJiYmFliqJKkDe7aac/H1pzh9ZvbLfPmtv2Cut5vGSZcgPAVs7pvfBJzsH1BV9/ZN35TkN5I8sarumTXuAHAAYHJycuVfsZJ0/rNnqzmnzxTTuzeMugyNoS6XRhwGtibZkmQdsAs41D8gyZclSW96e+9+Pz3sYiVJA9mzJamjgWeEq2o6yV7gZmANcF1VHU1ydW/9fuB7gR9JMg18AdhVVZ49kKQVZs+WpO4yqt43OTlZR44cGcm+JWkImro40J6t1WzttSe9NEKDzNmz/WY5SZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUmdgnCSHUnuSHIsyb55xj0ryf1Jvnd4JUqSFsKeLUndDAzCSdYA1wBXAtuAq5JsO8e41wM3D7tISVI39mxJ6q7LGeHtwLGqOl5VZ4GDwM45xv0o8GbgriHWJ0laGHu2JHXUJQhvBE70zU/1ln1Rko3Ai4D9wytNkrQI9mxJ6mhthzGZY1nNmv9V4Geq6v5kruG9O0r2AHsAJiYmutYoSerOnq2xdPH1pzh9ZvZLtZv1F5z7dSzNp0sQngI2981vAk7OGjMJHOw11CcCL0wyXVV/2D+oqg4ABwAmJycX92qXJM3Hnq2xdPpMMb17w6jLUGO6BOHDwNYkW4C/BXYBP9A/oKq2PDid5HeAP5ndUCVJK8KeLUkdDQzCVTWdZC8zv1m8Briuqo4mubq33mvMJGmVsGdLUnddzghTVTcBN81aNmczraqXLb0sSdJi2bMlqRu/WU6SJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqUqcgnGRHkjuSHEuyb471O5N8OMktSY4kuXz4pUqSurBnS1I3awcNSLIGuAa4ApgCDic5VFW39w17J3CoqirJpcANwFcvR8GSpHOzZ0tSd13OCG8HjlXV8ao6CxwEdvYPqKrPVVX1Zh8FFJKkUbBnS1JHXYLwRuBE3/xUb9mXSPKiJH8D/Cnwr4dTniRpgezZktTRwEsjgMyx7CFnD6rqLcBbkjwXeA3wgofcUbIH2AMwMTGxsEolSV3YszUSF19/itNnFv/hwvoL5nrpSsurSxCeAjb3zW8CTp5rcFW9J8lXJnliVd0za90B4ADA5OSkH8VJ0vDZszUSp88U07s3jLoMaUG6XBpxGNiaZEuSdcAu4FD/gCRflSS96WcC64BPD7tYSdJA9mxJ6mjgGeGqmk6yF7gZWANcV1VHk1zdW78f+B7gpUnuA74AfH/fL2JIklaIPVuSusuoet/k5GQdOXJkJPuWpCFo6oJGe7YGWXvtSS+N0Go2Z8/2m+UkSZLUJIOwJEmSmmQQliRJUpMMwpIkSWqSQViSJElNMghLkiSpSQZhSZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJnUKwkl2JLkjybEk++ZY/+IkH+7d/iLJ04dfqiSpC3u2JHUzMAgnWQNcA1wJbAOuSrJt1rCPA99UVZcCrwEODLtQSdJg9mxJ6q7LGeHtwLGqOl5VZ4GDwM7+AVX1F1V1ujf7fmDTcMuUJHVkz5akjtZ2GLMRONE3PwU8e57xrwD+bK4VSfYAewAmJiY6lihJWgB79ipx8fWnOH2mRl3Gill/QUZdgrRgXYLwXK/sOd/ZSb6ZmaZ6+Vzrq+oAvY/gJicn2+kOkrRy7NmrxOkzxfTuDaMuQ9I8ugThKWBz3/wm4OTsQUkuBa4FrqyqTw+nPEnSAtmzJamjLtcIHwa2JtmSZB2wCzjUPyDJBHAj8JKq+tjwy5QkdWTPlqSOBp4RrqrpJHuBm4E1wHVVdTTJ1b31+4GfB54A/EYSgOmqmly+siVJc7FnS1J3qRrNZV+Tk5N15MiRkexbkoagqd8Msmcv3NprT3qNsLR6zNmz/WY5SZIkNckgLEmSpCYZhCVJktQkg7AkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTeoUhJPsSHJHkmNJ9s2x/quT/GWSM0l+cvhlSpK6smdLUjdrBw1Isga4BrgCmAIOJzlUVbf3DfsM8GPAdy1LlZKkTuzZktRdlzPC24FjVXW8qs4CB4Gd/QOq6q6qOgzctww1SpK6s2dLUkddgvBG4ETf/FRvmSRp9bFnS1JHAy+NADLHslrMzpLsAfYATExMLOYuluzVH7z3i9O/cNlFI6lhocax5qVaymNe6vEa1fEe1/2O6rlq8X3R0XnVswEuvv4Up88s6iGM1PoL5noqJK0mXYLwFLC5b34TcHIxO6uqA8ABgMnJyZF0tdf89ee+OD0u/3iOY81LtZTHvNTjNarjPa77HdVz1eL7oqPzqmcDnD5TTO/eMKrdSzqPdbk04jCwNcmWJOuAXcCh5S1LkrRI9mxJ6mjgGeGqmk6yF7gZWANcV1VHk1zdW78/yZcBR4CLgAeS/ASwraruPecdS5KGzp4tSd11uTSCqroJuGnWsv1903/HzMdvkqQRs2dLUjd+s5wkSZKaZBCWJElSkwzCkiRJapJBWJIkSU0yCEuSJKlJBmFJkiQ1ySAsSZKkJhmEJUmS1CSDsCRJkppkEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUpE5BOMmOJHckOZZk3xzrk+TXe+s/nOSZwy9VktSFPVuSuhkYhJOsAa4BrgS2AVcl2TZr2JXA1t5tD/Dfh1ynJKmD865n3333l/4cF3ffDYcPj1/di7XUx7uU7Ud1rEf5HI/qeI3j8zRAlzPC24FjVXW8qs4CB4Gds8bsBH6vZrwfeFySpwy51qV705vmn1+NxrHmpVrKY17q8RrV8R7X/Y7quWrxfdHd+dWzL7lkZvqSS8bneX6w7iuuGK+6F2upj3cp24/qWI/yOR7V8RrH56mDVNX8A5LvBXZU1e7e/EuAZ1fV3r4xfwK8rqre15t/J/AzVXXkXPf7mMc8pi677LIhPISO7rsP3v9+3n3Rk+DMPwLw2C/cCxddBMnK1bEQVXDvvXz2kRd9cdGqr3mplvKYl3q8RnW8x3W/o3qu+re94EJ40ha+6dhfwXOeAw9/+OC6h+Rd73rXqnwTnm89mwce4H2/9bdc/sMb4WEPW/HnecH66v6icah7sZb6eJey/aiO9Sif41Edr3F8nmY5Z8+uqnlvwPcB1/bNvwT4b7PG/Clwed/8O4HL5rivPcCR3u22Qfse6g2eVfD3BfWbM/+UVm/+WStax/le8ygf81KP16iO97jud1TPVYvviwXc7Nkjvo1r3aN6vOPYC0b5HI9j313l74kuZ4S/HvgPVfVtvfmf7QXoX+wb85vAu6rqTb35O4DnVdWpDiF9xSU5UlWTo65jIcax5qVaymNe6vEa1fEe1/2O6rlq8X0xiD179RjXuhdrXPvIUozyOR7Hvrsa3xNdrhE+DGxNsiXJOmAXcGjWmEPAS3u/ifwc4LOrtaFK0nnOni1JHa0dNKCqppPsBW4G1gDXVdXRJFf31u8HbgJeCBwD/hF4+fKVLEk6F3u2JHU3MAgDVNVNzDTO/mX7+6YLeOVwS1tWB0ZdwCKMY81LtZTHvNTjNarjPa77HdVz1eL7YiB79qoxrnUv1rj2kaUY5XM8jn131b0nBl4jLEmSJJ2P/IplSZIkNanZIJzkR3tfQXo0yS+Nup4ukrym93WotyR5W5INo65pOSS5LsldSW5bxLabk/zvJB/tPbc/voBtH5HkA0lu7W376oXufymSfCLJR3rP7zn/nuuQ9/m03v4evN2b5CcWsP2/6R2r25K8KckjFrDtj/e2O9pln3O9LpI8Psnbk/zf3s/1Xfev8TNufdue3Wlbe/bC9jmynt3bvnPfHpuePeq/3zaKG/DNwDuAC3rzTxp1TR3rvqhv+seA/aOuaZke53OBZ7KIv1sKPAV4Zm/6McDHgG0dtw3w6N70w4G/Ap6zgo/7E8ATR3jc1wB/B1zScfxG4OPAI3vzNwAv67jt1wK3ARcy87sK7wC2LvR1AfwSsK83vQ94/aiOn7flvY1j37Znd9rWnr34/a9Yz+6NX1DfHpee3eoZ4R9h5luVzgBU1V0jrqeTqrq3b/ZRwHl5gXdVvQf4zCK3PVVVH+pN/wPwUWbe/F22rar6XG/24b3beXmMz+H5wP+rqk8uYJu1wCOTrGWmOZ7suN3XAO+vqn+sqmng3cCL5tvgHK+LncDv9qZ/F/iuroVr7Ixd37Znd9rWnr14K9mzYYF9e1x6dqtB+KnANyb5qyTvTvKsURfUVZLXJjkBvBj4+VHXs5ol+XLg65g5S9B1mzVJbgHuAt5eVZ23HYIC3pbkg0n2rOB+H7QL6PwF8FX1t8AvA3cCp5j5W7Rv67j5bcBzkzwhyYXM/CmvzQusF+DJ1fv7t72fT1rEfWg8jGXftmd3Z89esJXs2TCcvr3qevZ5G4STvKN3Hcvs205m/ke0HngO8FPADUnm/g7qFTagbqrq31XVZuD3gb2jrXb1SvJo4M3AT8w6KzOvqrq/qp4BbAK2J/na5apxDt9QVc8ErgRemeS5K7XjzHzxwncC/2sB26xn5n/3W4ANwKOS/GCXbavqo8DrgbcDbwVuBaYXWLbOM+PYt+3Zw2HPXpiV7tlw/vbtTn9HeBxV1QvOtS7JjwA31sxFKh9I8gDwRODularvXOare5Y/AP4U+IVlLGcsJXk4Mw3196vqxsXcR1X9fZJ3ATuY+V/wsquqk72fdyV5C7AdeM9K7JuZRv6hqvrUArZ5AfDxqrobIMmNwL8A/keXjavqjcAbe9v+Z2BqQRXP+FSSp1TVqSRPYeaskMbUOPZte/bS2bMXZcV7Ngylb6+6nn3enhEe4A+BbwFI8lRgHXDPSCvqIMnWvtnvBP5mVLWsVr0zRG8EPlpV/3WB216c5HG96Ucy0zRW5BgneVSSxzw4DXwrK9TMe65iAR+x9dwJPCfJhb3j/nxmru/rJMmTej8ngO9exP5h5quCf6g3/UPAHy3iPjQexq5v27MHs2cv2or3bBhK3151PbvJL9TofaRwHfAM4Czwk1X156OtarAkbwaeBjwAfBK4unfNz3klyZuA5zFztudTwC/0/hfaZdvLgfcCH2HmOAG8qma+aWvQtpcyc/H+Gmb+k3hDVf3HBT+ARUjyFcBberNrgT+oqteu0L4vBE4AX1FVn13gtq8Gvp+Zj8f+Gtj94C8zddj2vcATgPuAf1tV7xww/iGvC2bC0Q3ABDNN/vuqalG/tKPVbRz7tj2707b27IXveyQ9u7d95749Lj27ySAsSZIktXpphCRJkhpnEJYkSVKTDMKSJElqkkFYkiRJTTIIS5IkqUkGYUmSJDXJICxJkqQmGYQlSZLUpP8PC8pPkGAB+WoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# nbi:hide_in\n",
    "# library\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "plt.rcParams[\"figure.figsize\"] = (12, 5)\n",
    "\n",
    "K=10\n",
    "range_x=np.array([-3, -6, 5, 9, 2, 1, 3, 8, 7, 10])\n",
    "pmf_values=np.ones(range_x.size)/range_x.size\n",
    "\n",
    "# Sort x and y according to order in x\n",
    "sortargs = range_x.argsort()\n",
    "range_x = range_x[sortargs]\n",
    "pmf_values = pmf_values[sortargs]\n",
    "\n",
    "#cdf values using cumsum function with padding\n",
    "cdf_values = np.cumsum(pmf_values)\n",
    "\n",
    "def padding(cdf_values, range_x):\n",
    "    edge = (range_x[0]-2, range_x[-1]+2)\n",
    "    range_x_padded = np.pad(range_x, (1,1), 'constant', constant_values=edge)\n",
    "    cdf_values_padded = np.pad(cdf_values, (1,1), 'constant', constant_values=(0, 1))\n",
    "    return cdf_values_padded, range_x_padded, edge\n",
    "\n",
    "cdf_values_padded, range_x_padded, edge = padding(cdf_values, range_x)\n",
    "# setting up the figure\n",
    "fig, [ax1, ax2] = plt.subplots(1,2, num=1, clear=True)\n",
    "\n",
    "ax1.set_ylim(0,1.1)\n",
    "ax1.axhline(y=0, color='k')\n",
    "ax1.set_xlim(edge)\n",
    "ax1.set_xticks(range_x)\n",
    "ax1.set_yticks(cdf_values_padded)\n",
    "ax1.spines[\"top\"].set_visible(False)  \n",
    "ax1.spines[\"right\"].set_visible(False)\n",
    "ax1.spines[\"bottom\"].set_visible(False)\n",
    "\n",
    "ax2.set_ylim(0,1.1) \n",
    "ax2.axhline(y=0, color='k')\n",
    "ax2.set_xlim(edge)\n",
    "ax2.set_xticks(range_x)\n",
    "ax2.set_yticks(cdf_values_padded)\n",
    "ax2.spines[\"top\"].set_visible(False)  \n",
    "ax2.spines[\"right\"].set_visible(False)\n",
    "ax2.spines[\"bottom\"].set_visible(False)\n",
    "\n",
    "# plot line grapghs\n",
    "ax1.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20)\n",
    "markers,stems,base = ax1.stem(range_x, pmf_values, markerfmt=' ',\n",
    "                              linefmt='#039be5', basefmt=\"black\", use_line_collection=True)\n",
    "stems.set_linewidth(3)\n",
    "ax1.set_title(\"pmf\")\n",
    "\n",
    "# plot cdf using step function\n",
    "ax2.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20)\n",
    "ax2.step(range_x_padded, cdf_values_padded, linewidth=1, where='post', color='#039be5',)\n",
    "ax2.set_title(\"cdf\")\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Relative frequency histogram for pmf\n",
    "\n",
    "Suppose our data   consists of values of the random variable on a sequence of trial outcomes. As mentioned earlier, probability of an event represents how frequent the experiment outcome terminates in the event, in a large number of repetitive trials. Hence,  the pmf at $x\\in\\text{Range}(X)$ can be empirically estimated using the relative frequency \n",
    "\n",
    " $$f_{\\text{emp}}(x)= \\frac{\\text{number of elements in data = }x}{ \\text{ size of data}}.$$\n",
    " \n",
    "Resulting relative frequency histogram will approximate the pmf histogram.\n",
    "\n",
    "In the example below, we consider the following experiment: a dice is tossed twice and the random variable is the maximum of the two tosses: \n",
    "\n",
    "$$S=\\{(i,j):\\; 1\\leq i\\leq n, 1\\leq j\\leq 6\\}$$\n",
    "\n",
    "and for any $s=(i,j)$,\n",
    "      \n",
    " $$X(i,j)=\\max\\{i,j\\}.$$\n",
    " \n",
    " It can be seen that in this case, $\\text{range}(X)=\\{1,\\dots, n\\}$ and, for any $x\\in \\text{range}(X)$,\n",
    " \n",
    " $$f(x)=\\frac{2x-1}{n^2}.$$\n",
    "\n",
    "The data consist of values of $X$ on 1000 random pairs of tosses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAEzCAYAAAD+XEDdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de3hU1bk/8O+bgKUFqhURclGhFEq5JDEZI/5QQUSBkoC1KGCtSg8EfQ76a1FbrVqS9vm1whFOaw8tCRy1p1rirXAIxUttBUsvliRFCioWKWgyEVIUJFUKSd7fH5lMZpLZe3aA7LWT9f08zzxk9mX2u1529puZvWYtUVUQERFRMKWYDoCIiIicsVATEREFGAs1ERFRgLFQExERBRgLNRERUYCxUBMREQWYp0ItIlNFZLeI7BGRe1y2u0hEmkRkVmf3JSIioo6SFmoRSQWwEsA0AKMAzBWRUQ7bLQXwQmf3JSIiosS8vKPOB7BHVfeq6nEA5QBmJtjudgDPAjh4EvsSERFRAl4KdQaAd2Oe10SWRYlIBoAvAVjV2X2JiIjIWS8P20iCZe3HHf0hgG+papNI3OZe9m3ZUKQIQBEAjBo1Km/Xrl0eQiMiIuoREtVLAN4KdQ2A82KeZwIIt9smBKA8UqTPAfBFEWn0uC8AQFXLAJQBQCgU4gDkRERE8FaotwEYLiJDAdQCmAPghtgNVHVo688i8hiAjaq6XkR6JduXiIiInCUt1KraKCKL0NKbOxXAI6q6S0Rujaxvf1866b6nJ3QiIqKeT4I4zWUoFNLKykrTYRAREfnllO5RExH56sSJE6ipqcGxY8dMh0J0WvXp0weZmZno3bu3531YqIkocGpqatC/f38MGTIE7b5JQtRtqSoOHTqEmpoaDB06NPkOERzrm4gC59ixYxgwYACLNPUoIoIBAwZ0+pMiFmoiCiQWaeqJTua8ZqEmImrn0KFDyMnJQU5ODgYPHoyMjIzo8+PHj5+246xZswZf//rXE66bMmUKjh496rjvihUreA/fErxHTUSBN6z8APY3NJ2217ugXyrenjPIcf2AAQOwfft2AEBxcTH69euHu+66K24bVYWqIiWla97vvPDCC67rV6xYga997Wvo06fPKR+rsbERvXqxHAQV/2eIKPD2NzShcX76aXu9XmsSDpCY1J49e3DNNdfg0ksvxauvvor169cjOzsbhw8fBgCUl5fjpZdewpo1a3DgwAHcdttteOedd5CSkoKHH34Y48aN6/CaNTU1mDJlCvbu3YtZs2bhBz/4AQAgMzMTO3fuRGpqKq6//nqEw2E0NTWhuLgY7777Lg4ePIjLLrsMgwYNwksvvYTHH38cS5cuhapixowZ+P73vw8AKC0txfLly5Geno7Pfe5z6NevH374wx/ixhtvxKBBg1BdXY2LLroI1157Lb7xjW/g2LFj+NSnPoXHHnsMw4cPx5o1a/D888/j6NGjHWIkf7BQExF1wuuvv45HH30Uq1atQmNjo+N2d9xxB775zW9i3Lhx2LdvHwoKCrBz584O27322muorq5Gr169MGLECNx+++1IT2/7o2TTpk0YMmQInnvuOQDAkSNHcOaZZ2L58uX43e9+h7POOgs1NTW4//77UVlZiTPPPBOTJ0/Gxo0bkZ2djQcffBDV1dXo27cvJk6ciPz8/Ohrv/322/jNb36DlJQUHDlyBFu3bkVqaiqef/553H///XjyySc9xUhdi4WaiKgThg0bhosuuijpdi+99BJ2794dff7BBx/g448/xic/+cm47SZPnoz+/fsDAEaOHIl33nknrghmZWXhnnvuwT333IPCwkKMHz++w7FeffVVTJo0Ceeccw4A4IYbbsArr7yCY8eOYdKkSfjMZz4DAJg1axbeeeed6H7XXXdd9KP7w4cP46abbsLbb7/d4fWTxUhdi53JiIg6oW/fvtGfU1JSEDu6Y2znLlXFn//8Z2zfvh3bt29HbW1thyINAJ/4xCeiP6empnZ4l/6FL3wBlZWVGD16NO6+++7oR9qxnEaYTDbyZGxb7rvvPkyZMgU7d+7E+vXr49qSLEbqWizUREQnKSUlBZ/5zGfwt7/9Dc3NzVi3bl103eTJk7Fy5cro89bOaZ1VW1uLfv364atf/SoWL16M6upqAED//v2jvcLHjRuHl19+GYcOHUJjYyPKy8sxYcIEXHzxxXj55Zdx+PBhnDhxAr/85S8dj3PkyBFkZGQAAB577LGTipW6Bgs1EdEpWLp0KaZOnYorr7wSmZmZ0eUrV67E73//e2RlZWHUqFFYvXr1Sb3+a6+9hosuugg5OTlYtmwZvv3tbwMAioqKMHnyZEyePBmZmZn47ne/i4kTJyInJwfjxo3D9OnTcf755+Puu+9Gfn4+rr76aowePRpnnnlmwuN861vfwt13353wo3Uyi5NyEFHgvPHGG/jCF74Qfe7317N6koaGBvTr1w8nTpzAzJkzcdttt6GwsNB0WFZrf35HcFIOIuq+bCmqXeGBBx7A5s2bcezYMUydOhUFBQWmQ6JOYqEmIurB/vM//9N0CHSKeI+aiIgowFioiYiIAoyFmoiIKMBYqImIiAKMhZqIKIHU1FTk5ORgzJgxKCwsjE684aZfv36u6w8fPoyf/OQn0efhcBizZs065VgB4O67746OXtZT3HLLLXjmmWc6LE+Wt/Z57u5YqImoWyiqKIKUSPQRPhpGxe6KuGVlVWUAELescG3Ld4YL1xbGLU/mk5/8JLZv346dO3fi7LPPjhtl7GS1LyDp6ekJC9HJKC0tRXV1Nf7jP/4jbnlPHO4zWd5Od6Fuajp93+E/Ka1zqgbpkZeXp0Rkr9dffz3ueW5pru8x9O3bN/rzT3/6U73tttuiz5ctW6ahUEjHjh2r3/nOdzrsc/ToUZ00aZJeeOGFOmbMGF2/fr2qqs6ePVv79Omj2dnZetddd+nf//53HT16tKqq5ufn686dO6OvNWHCBK2srNSGhgadN2+ehkIhzcnJib5WrMLCQk1JSdHs7GwtLy/Xm2++Wb/xjW/oxIkTdfHixY6v8dFHH+ns2bN17Nixev3112t+fr5u27atQ/uffvppvfnmm1VV9eDBg3rttddqKBTSUCikW7duVVXVJUuW6Lx583TChAk6dOhQ/dGPfhTd/2c/+5mOHTtWs7Ky9MYbb9QPP/xQhwwZosePH1dV1SNHjugFF1wQfd7q5ptv1ttvv10vueQSHTp0qD799NOqqnF527lzp1500UWanZ2tY8eO1bfeeqtDnpubm/Wuu+7S0aNH65gxY7S8vFxVVZuamvS2227TUaNG6fTp03XatGnRY1xwwQVaUlKi48eP17Vr12pZWZmGQiHNysrSa6+9Vv/5z39GY7z11lt14sSJOnToUN28ebPOmzdPR44cGc1Ze+3P7wjHmmi8KCd6sFAT2a39hQzFOK2vX1pZmnSb1kLV2Nios2bN0ueee05VVV944QVdsGCBNjc3a1NTk06fPl23bNkSt8+JEyf0yJEjqqpaX1+vw4YN0+bm5rgCoxpfcFasWBEt+uFwWIcPH66qqvfee6/+/Oc/V1XVDz74QIcPH64NDQ2O8aq2FI/p06drY2Oj62ssX75c582bp6qqr732mqampiYt1HPnztXf/e53qqq6f/9+HTlypKq2FOpLLrlEjx07pvX19Xr22Wfr8ePHdefOnTpixAitr69XVdVDhw6pquott9yi69ata/n/KC3VxYsXd2jTzTffrLNmzdKmpibdtWuXDhs2rEPeFi1apI8//riqqv7rX//Sjz76qEOen3nmGZ08ebI2Njbqe++9p+edd56Gw2F9+umnddq0adrU1KR1dXV61llnxRXqpUuXRl/jH//4R/Tn++67Tx9++OFojLNnz9bm5mZdv3699u/fX3fs2KFNTU2am5urf/nLXzq0q7OFmh99E5F1Fm5cmHSbjz/+GDk5ORgwYADef/99XHXVVQCAF198ES+++CIuvPBC5Obm4s0338Tf/va3uH1VFd/+9reRlZWFyZMno7a2FgcOHHA93vXXX4+nn34aAPDUU0/huuuuix7vwQcfRE5ODiZOnIhjx47FTVXp5LrrrkNqaqrra7zyyiu48cYbAbRMp5mVlZX0dV966SUsWrQIOTk5mDFjBj788MPo5CDTp0/HJz7xCZxzzjk499xzceDAAfz2t7/FrFmzolNwnn322QCA+fPn49FHHwUAPProo5g3b17C411zzTVISUnBqFGjEubwkksuwfe//30sXboU+/fvTzhD2datWzF37lykpqZi0KBBmDBhArZt24atW7dGp/ocPHgwrrjiirj9Zs+eHf15586duOyyyzB27Fg88cQT2LVrV3RdYWEhRARjx47FoEGDMHbsWKSkpGD06NHYt29f0pwmw5HJiCjw0vql+X7M1nvUR44cQUFBAVauXIk77rgDqop7770XCxc6F/snnngC9fX1qKqqQu/evTFkyJC4aSMTycjIwIABA7Bjxw48+eSTKC0tBdBS9J999ll8/vOf71T8sVNYur2GSOL79bHLY2Nvbm7GH//4R89TdqpqwmOMHz8e+/btw5YtW9DU1IQxY8YkjCP2NTXB3BQ33HADLr74YvzqV7/ClClTsGbNGnz2s5+N2ybRfm7LW8Xm8JZbbsH69euRnZ2Nxx57DJs3b+4QY0pKSly8KSkpp6WPgKd31CIyVUR2i8geEbknwfqZIrJDRLaLSKWIXBqzbp+I/LV13SlHTETWCd8ZNnbsM888Ew8//DAeeughnDhxAlOmTMEjjzyChoYGAC3TUB48eDBunyNHjuDcc89F79698fLLL2P//v0A4qemTGTOnDlYtmwZjhw5grFjxwIApkyZgh//+MfRovKXv/yl021weo3LL78cTzzxBICWd4w7duyI7jNo0CC88cYbHabvvPrqq/Ff//Vf0efJpu+88sor8dRTT+HQoUMAgPfffz+67qabbsLcuXMd3017sXfvXnz2s5/FHXfcgRkzZmDHjh0d8nz55ZfjySefRFNTE+rr6/HKK68gPz8fl156KZ599lk0NzfjwIEDccW3vaNHjyItLQ0nTpyI5swvSQu1iKQCWAlgGoBRAOaKyKh2m/0GQLaq5gD4GoA17dZfoao5qho6DTETkWWKNxef1tfbMGdDp7a/8MILkZ2djfLyclx99dW44YYbcMkll2Ds2LGYNWtWh+L7la98BZWVlQiFQnjiiScwcuRIAMCAAQMwfvx4jBkzJuHXqGbNmoXy8nJcf/310WUPPPAATpw4gaysLIwZMwYPPPBAp9vr9Bq33XYbGhoakJWVhWXLliE/Pz+6z4MPPoiCggJMmjQJaWltn2g8/PDDqKysjE7fuWrVKtdjjx49Gvfddx8mTJiA7OxsLF68OC5PH3zwAebOndvpNrV68sknMWbMGOTk5ODNN9/ETTfd1CHPX/rSl5CVlYXs7GxMmjQJy5Ytw+DBg/HlL38ZmZmZGDNmDBYuXIiLL77YcRrQ733ve7j44otx1VVXRf8//ZJ0mksRuQRAsapOiTy/FwBU9Qcu2z+iql+IPN8HIKSq//AaFKe5JLJb+2kA23+dqnJBy/UhtLrtb/8lE5ageGIx0peno66hDgCQm5aLqqIqFFUUYXV123zQtYtrkd4/vSub0C1NnDgRDz30EEIhf95TPfPMM/jf//1f/PznP/fleIm0TgN66NAh5Ofn4/e//z0GDx7cpcfsimkuMwC8G/O8BsDFHY4g8iUAPwBwLoDpMasUwIsiogBKVbXMwzGJiKJ0icM9xgTLE31MXlZYhrJCXnqC5Pbbb8dzzz2HTZs2GY2joKAAhw8fxvHjx/HAAw90eZE+GV4KdaIq3+G3Q1XXAVgnIpcD+B6AyZFV41U1LCLnAvi1iLypqq90OIhIEYAiADj//PO9xk9ERKeJ2z3a0+3HP/6xb8dy42ebT5aXzmQ1AM6LeZ4JwLFnR6QIDxORcyLPw5F/DwJYByDfYb8yVQ2pamjgwIEewyciIurZvBTqbQCGi8hQETkDwBwAcT0xRORzEul/LyK5AM4AcEhE+opI/8jyvgCuBrDzdDaAiHqmZP1niLqjkzmvk370raqNIrIIwAsAUtHSUWyXiNwaWb8KwJcB3CQiJwB8DGC2qqqIDELLx+Gtx/qFqj7f6SiJyCp9+vTBoUOHMGDAAMfv+RJ1N6qKQ4cOoU+fPp3aL2mvbxPY65vIbidOnEBNTU3SQUKIups+ffogMzMTvXv3br/qlHp9ExH5qnfv3hg6dKjpMIgCgWN9ExERBRgLNRERUYCxUBMREQUYCzUREVGAsVATEREFGAs1ERFRgLFQExERBRgLNRERUYCxUBMREQUYCzUREVGAsVATEREFGAs1ERFRgLFQExERBRgLNRERUYCxUBMREQUYCzUREVGA9TIdABERdS/ho2FkrMiIPl+QuwBlhWXIK8tDdV01ACCtXxrCd4ZRvLkYJVtKottWLqgEAIRWh6LLlkxYguKJxUhfno66hjoAQG5aLqqKqlBUUYTV1auj29YurkVVuAozymdEl5UWlKIorwhSItFlBSMKUDG3AoVrC7HxrY3R5bpEUVZVhoUbF0aXbZizAXnpea5tuuDMC7Dv6/tOOmenQlTVyIHdhEIhraysNB0GERElULG7AoWfLzQdhq/Sl6cjfGe4Kw8hTiv40TcREXVK7LtZW3RxkXbFQk1ERJRE8eZiY8dmoSYiIkoi9j6731ioiYioU0oLSk2HYBUWaiIi6pSivCLTIViFhZqIiDol9mtQtmj9WpkJnr5HLSJTAfwIQCqANar6YLv1MwF8D0AzgEYAX1fVrV72JSKi4BtWfgD7G5qiz3utMdcL2gQ9Xg85o63NF/RLxdtzBvly7KSFWkRSAawEcBWAGgDbRGSDqr4es9lvAGxQVRWRLABPARjpcV8iIgq4/Q1NaJyfDgCQEkR/toWUZECXtI074ucfKl4++s4HsEdV96rqcQDlAGbGbqCqDdo2ckpfAOp1XyIi6l4KRhSYDsEqXgp1BoB3Y57XRJbFEZEvicibAH4F4Gud2ZeIiLqPirkVpkOwipdCnajXQIdxR1V1naqOBHANWu5Xe94XAESkSEQqRaSyvr7eQ1hERGRC4Vq7hg8FWsYjN8VLoa4BcF7M80wAjh/Oq+orAIaJyDmd2VdVy1Q1pKqhgQMHegiLiIhMiJ3kwhbFE4uNHdtLod4GYLiIDBWRMwDMAbAhdgMR+ZyISOTnXABnADjkZV8iIqKgS19urvNc0l7fqtooIosAvICWr1g9oqq7ROTWyPpVAL4M4CYROQHgYwCzI53LEu7bRW0hIiLqEq3Tb5rg6XvUqroJwKZ2y1bF/LwUwFKv+xIRUfcV+zUl6nocmYyIiDqlrKrMdAi+y03LNXZsFmoiIuqUhRsXmg7Bd1VFVcaOzUJNRESURFGFuYlIWKiJiIiSWF292tixWaiJiKhTNszht2z9xEJNRESdkpeeZzoEq7BQExFRp2SssG/KhtrFtcaOzUJNRESURFWYvb6JiIgCa0b5DGPHZqEmIqJOWZC7wHQIVmGhJiKiTikrtG9kMpNYqImIqFPyyuzr9V1aUGrs2CzURETUKdV11aZD8F1RHkcmIyIiCiwpEWPHZqEmIqJOSeuXZjoEq7BQExFRp4TvDJsOwSos1ERE1CnFm4tNh+C7ghEFxo7NQk1ERJ1SsqXEdAi+q5hbYezYLNRERERJFK4tNHZsFmoiIqIkNr610dixWaiJiKhTKhdUmg7BKizUREREAcZCTUREnRJaHTIdgu90iRo7Ngs1ERFREmVV5iYiYaEmIiJKYuHGhcaO7alQi8hUEdktIntE5J4E678iIjsijz+ISHbMun0i8lcR2S4i7IFARNTNLZmwxHQIVumVbAMRSQWwEsBVAGoAbBORDar6esxmfwcwQVU/EJFpAMoAXByz/gpV/cdpjJuIiAwpnlhsOgSreHlHnQ9gj6ruVdXjAMoBzIzdQFX/oKofRJ7+CUDm6Q2TiIiCIn15uukQfLdhzgZjx/ZSqDMAvBvzvCayzMm/AXgu5rkCeFFEqkTE3ISeRER0WtQ11JkOwXd56XnGjp30o28AiSbhTNhPXUSuQEuhvjRm8XhVDYvIuQB+LSJvquorCfYtAlAEAOeff76HsIiIus6w8gPY39BkOgwKiIwVGca+ouXlHXUNgPNinmcC6DDHmYhkAVgDYKaqHmpdrqrhyL8HAaxDy0fpHahqmaqGVDU0cOBA7y0gIuoC+xua0Dg/Pe7xkws3oqk2I/pYd1kV3pmNuGVfG1SMxvnpyG4ujC4790gIjfPTcf/nyuK2ffWLdXj1i3Vxy+7/XBka56fj3COh6LLs5kI0zk/H1wYVx237zmxg3WVVcct+cuFGNM5Pj1s2re9CNM5Px7S+C+OWJ2tTZsP/ibY9Vm5arqH/FTuJqvtfCCLSC8BbAK4EUAtgG4AbVHVXzDbnA/gtgJtU9Q8xy/sCSFHVo5Gffw3gu6r6vNsxQ6GQVlaygzgRmdNrTTiuQEmJGB30woTYNrfPh23a//93QT4SfXoNwMM7alVtBLAIwAsA3gDwlKruEpFbReTWyGbfATAAwE/afQ1rEICtIvIagD8D+FWyIk1ERMGwIHdBwuVFFfZ1N3LKhR+SvqM2ge+oicg0vqOOF5sP23MBBOwdNRERAQUjCkyH4Lu8MnM9nYPGZC5YqImIPKiYW2E6BN9V11WbDiEwTOaChZqIyIPCtYWmQwiM2sW1pkOwCgs1EZEHG9/aaDoE36X1S0u4vCpc5XMk5jnlwg8s1ERElFD4zg5DZgAAZpTP8DkS85xy4QcWaiIiSqh4c7HpEALDZC5YqImIPLDx60glW0pMhxAYJnPBQk1E5EFZVZnpEAKjtKDUdAhWYaEmIvJg4caFpkMIjKI8+0YmM4mFmoiIEqpckHiESClxHESrx3LKhR9YqImIiAKMhZqIyIMNczaYDsF3odUh0yEEhslcsFATEXmQl85xr1vZOO65SSzUREQeZKzIMB1CYNg47rlJLNRERJTQkglLEi63cdxzp1z4gYWaiIgSKp5YnHC5jeOeO+XCDyzUREQeLMhdYDoE36UvTzcdQmCYzAULNRGRB2WF9o1MVtdQZzqEwDCZCxZqIiIP8srY67uVjeOem8RCTUTkQXVdtekQfJeblptwuY3jnjvlwg8s1ERElFBVUVXC5TaOe+6UCz+wUBMReZDWL810CL4rquDkG61M5oKFmojIg/CdYdMh+G519WrTIQSGyVywUBMReVC8udh0CIFh47jnJrFQExF5ULKlxHQIgcFxz/3lqVCLyFQR2S0ie0TkngTrvyIiOyKPP4hIttd9iYgomGoX1yZcbuO450658EPSQi0iqQBWApgGYBSAuSIyqt1mfwcwQVWzAHwPQFkn9iUiogCqCpvr6Rw0JnPh5R11PoA9qrpXVY8DKAcwM3YDVf2Dqn4QefonAJle9yUi6g4qF1SaDsF3M8pnmA4hMEzmwkuhzgDwbszzmsgyJ/8G4LmT3JeIiALOxnHPTfJSqCXBsoTjx4nIFWgp1N86iX2LRKRSRCrr6+s9hEVE5J/Q6pDpEALDxnHPTfJSqGsAnBfzPBNAhy8UikgWgDUAZqrqoc7sCwCqWqaqIVUNDRw40EvsRETUhUoLShMut3Hcc6dc+MFLod4GYLiIDBWRMwDMARD3JToROR/ALwF8VVXf6sy+REQUTEV5iUfjsnHcc6dc+CFpoVbVRgCLALwA4A0AT6nqLhG5VURujWz2HQADAPxERLaLSKXbvl3QDiKiLrVkwhLTIfhOShLdvbSTyVz08rKRqm4CsKndslUxP88HMN/rvkQUTMPKD2B/Q5PpMAKpeGKx6RACw8Zxz03yVKiJyA77G5rQOD/ddBiB0GtNfHea9OXpVo73nQjz4C8OIUpE5EFdQ53pEHxXMKIg4XIbxz13yoUfWKiJiCihirkVCZfbOO65Uy78wEJNRORBblqu6RB8V7i20HQIgWEyFyzUREQeVBXZN+71xrc2mg4hMEzmgoWaiMiDogpz36MNGhvHPTeJvb6JKKHY740WjChAxdwKFK4tjHtnoUsUZVVlWLhxYXTZhjkbkJeeFzcV4oLcBSgrLENeWV50sIy0fmkI3xlG8ebiuHuerUUgdsjOJROWoHhiMdKXp0c7deWm5aKqqApFFUVYXb06um3t4lpUhaviJlEoLShFUV6R5zYlsrp6NYfOJCNENfFJaVIoFNLKSv7FRuS3XmvC/HpWRPtcSIk4FnEbxObD9lwAXfK74jiiCj/6JqKEyqrse/doY5vdMB9tTOaChZqIEor9ONsWbm2uXVzrYyTBYOM54MRkLlioiYg8qArb1+vbiY3jnpvEQk1E5EFs5zTbcdxzf7FQE1FCG+bYNyOtjW1245SP9OX2dTg0eW6wUBNRQnnpeaZD8J2NbXbjlA8bxz03eW6wUBNRQrHfg7aFW5tLC0p9jCQYbDwHnJjMBQs1EZEHRXkcmayVjeOem8RCTUTkQeyoZrazcdxzk1ioiSihBbkLTIfgOxvb7MYpHzaOe27y3GChJqKEbBzX2sY2u3HKR+zY6rYweW6wUBNRQnll9vWAdmtzwYgCHyMJBhvPAScmc8FCTUQJtc5yZRO3NlfMrfAxkmCw8RxwYjIXLNRERB4Uri00HUJg2DjuuUks1ESUUFq/NNMh+M6tzbFzVtvCKR82jntu8veBhZqIEgrfGTYdgu9sbLMbp3zYOO65yXODhZqIEireXGw6BN/Z2GY3zEcbk7nwVKhFZKqI7BaRPSJyT4L1I0XkjyLyLxG5q926fSLyVxHZLiKVpytwIupaJVtKTIfgO7c26xL1MZJgsPEccGIyF0kLtYikAlgJYBqAUQDmisiodpu9D+AOAA85vMwVqpqjqqFTCZaIyJSyKn7HupWN456b5OUddT6APaq6V1WPAygHMDN2A1U9qKrbAJzoghiJiIxbuHGh6RACg+Oe+8tLoc4A8G7M85rIMq8UwIsiUiUi/N8l6iYqF9h3p8rGNrtxyoeN456bPDd6edgm0f9IZ27WjCuasBMAABRSSURBVFfVsIicC+DXIvKmqr7S4SAtRbwIAM4///xOvDwREVHP5eUddQ2A82KeZwLw3E9dVcORfw8CWIeWj9ITbVemqiFVDQ0cONDryxNRFwmttq9LiVubN8zZ4GMkwWDjOeDEZC68FOptAIaLyFAROQPAHACezlgR6Ssi/Vt/BnA1gJ0nGywRkSl56Rz3upWN456blPSjb1VtFJFFAF4AkArgEVXdJSK3RtavEpHBACoBfBpAs4h8HS09xM8BsE5EWo/1C1V9vmuaQkTUdTJWZFj5Fa1EbBz33CQv96ihqpsAbGq3bFXMz++h5SPx9j4EkH0qARKRGUsmLDEdgu9sbLMbp3wUri20rlibPDc4MhkRJVQ8sdh0CL6zsc1unPJh47jnJs8NFmoiSih9ebrpEHzn1uYFuQt8jCQYbDwHnJjMhaePvol6qmHlB7C/ocl0GIFU11BnOgTfubW5rNC+kclsPAecmMwFCzVZbX9DExrn811Dq15rOHuUk7yyPFQV2Te9YyLsVOcvfvRNFGPID4dEZ8lJX54OKRFIiSCvrOWrOUUVRdFlUiIIHw2jYndF3LLWMaFjlxWuLQTQ0gkndjnQMoZ07LKK3RUIHw3HLSuqaBnUL68sL7qs9aO44s3FcdtWhatQFa6KW+bWJie5abmnNbfdgVubq+uqfYwkGJzyYeO45yZ/H0Q1eH8ZhUIhrazkUH7U9XqtCce9o5YSserdQlFFUdxHuu3zYTPbz432YvNhey6ALvldcRyXle+oiSy2unq147rWd/E2cWtzWr80HyMJBhvPAScmc8FCTRSDkzK0cSviPZVbm8N32nf/3sZzwInJXLBQExF50Hqfn+wc99wkFmqiGLZNQlC7uNZ0CN1GyZYS0yEEBsc99xcLNZHFqsLOXzeysYjb2GY3TvnIWJHhcyTmmTw3WKiJLDajfIbjOrci3lPZ2GY3zEcbk7lgoSaKwUkZ2rgV8Z7Krc02djS08RxwYjIXLNREMTgpA1FyNo57bhILNVEM2yYhKC0oNR1Ct2FbR0M3No57bhILNVEM2yYhKMpzHsTBxiJuY5vdOOUj2fCzPZHJc4OFmshireONJ+JWxHsqG9vsxikfNo57bvLcYKEmimHjRBRO3Ip4T+XWZhs7Gtp4DjgxmQsWaqIYnMaQnLCjYRsbxz03iYWaKIZtkxAUjCgwHUK3YVtHQzc2jntuEgs1UQzbJiGomFvhuM7GIu7WZts6GgLO+bBx3HOTvw8s1EQWK1xb6LjOrYj3VDa22Y1TPmwc99zkucFCTWSxjW9tdFznVsR7Krc229jR0MZzwInJXLBQE8XgpAxt3Ip4T+XWZhs7Gtp4DjgxmQsWaqIYnISAnNjW0dCNjeOem+SpUIvIVBHZLSJ7ROSeBOtHisgfReRfInJXZ/YlChLbJiHQJWo6hG7Dto6GFBxJC7WIpAJYCWAagFEA5orIqHabvQ/gDgAPncS+RGRIWZXzmM02FnEb2+zGKR82jntu8tzw8o46H8AeVd2rqscBlAOYGbuBqh5U1W0ATnR2XyIyZ+HGhY7r3Ip4T2Vjm90wH21M5sJLoc4A8G7M85rIMi9OZV8i33FShjZuRbyncmuzjR0NbTwHnJjMhZdCnWiAU6+fAXjeV0SKRKRSRCrr6+s9vjzR6cVJGcgJOxq2sXHcc5O8FOoaAOfFPM8E4HX8OM/7qmqZqoZUNTRw4ECPL090etk2CcGGORtMh9Bt2NbR0A3HPfeXl0K9DcBwERkqImcAmAPA62/3qexLRF0sL915XmEbi7iNbXbjlA8bxz03eW4kLdSq2ghgEYAXALwB4ClV3SUit4rIrQAgIoNFpAbAYgD3i0iNiHzaad+uagwRdU7GCucuI25FvKeysc1unPJh47jnJs8NT9+jVtVNqjpCVYep6v+LLFulqqsiP7+nqpmq+mlVPSvy84dO+xIFlY0TUThxK+I9lVubbexoaOM54MRkLjgyGVEMTspATtjRsI2N456bxEJNFMO2SQgW5C4wHUK3YVtHQzc2jntuEgs1UQzbJiEoK3QexMHGIm5jm9045cPGcc9Nnhss1EQWyytz7iDjVsR7Khvb7MYpHzaOe27y3GChJrJYdV214zq3It5TubXZxo6GNp4DTkzmgoWaKAYnZWjjVsR7Krc229jR0MZzwInJXPQydmQyZlj5AexvaAIAaNN7aH6v7S9F+dRXkPKZZWg6OBU48deWhSmDkJpWjeYPl0OProhumzLwOQBAc/20tv37L0bKp+9EU10u0HygZWHvsUg993k0f/BN6EdPtO0/uAo4vgPN789r2/+spUjpeyOaamO+CtFnMlIH/AxNh24Gjr0UXZyaUYvmfz4OPfytttc8+1HgjCz3NqVmInXwqwlzU1ZVZlXv3rR+aaZD6DYK1xZaWawTsXHcc5NYqC20v6EJjfNbRhaq2F2Fws8nehe5I8Gy5ZFHe4n2fy/Bsscjj1i5AG7x+Jq/TrDsm5GHl/1b2pS+PB3hSPt7rYkf0XbhxoVWFerwnc6jAdtYxN3abFtHQ8A5H1XhKqR/3q7RyUz+PvCjb8vZOH6xW3GyTfHmYsd1NubJxja7ccoHrxv+YqEm67gVJ9uUbClxXGdjnmxssxvmo43JXLBQk3XcihMnZWjjlqeeyq3NNnY0tPEccGIyFyzUlrNx/GI3nJSBnJRV8TvWrXjd8BcLteVs6jjlhW2TEFQuqDQdQrexcONC0yEEBq8b/mKhtpyN4xezOHljY55sbLMbp3zwuuEvFmoii4VWh0yHQERJsFCTddyKEydlaGNjEXdrs40dDW08B5yYzAULteVsHL/YDSdlICfsaNiG1w1/sVBbjkMixrNtEoIlE5aYDqHbsK2joRteN/zFQm25wrWFpkPwnVtxsm0SguKJxY7rbCziNrbZjVM+eN3wFwu15Wwcv9itONkmfbnzeM025snGNrtxygevG/5ioSbruBUn2yaiqGuoc1znlqeeyq3NNnY0tPEccGIyFyzUZB234sRJGdq45amncmuzjR0NbTwHnJjMBQu15Wwcv9iNbZMQ5Kblmg6h27Cto6EbXjf8xUJtORvHL3YrTrZNQlBVVOW4zsYi7tZm2zoaAs754HXDXyzUlrNx/GK34mSbogrnMZttzJONbXbjlA9eN/zlqVCLyFQR2S0ie0TkngTrRUQejqzfISK5Mev2ichfRWS7iHAgXTLOrTjZZnX1asd1NubJrc22dTQE7DwHnJjMRdJCLSKpAFYCmAZgFIC5IjKq3WbTAAyPPIoA/LTd+itUNUdVOR4dGedWnDgpQxu3PPVUbm22saOhjeeAE5O58PKOOh/AHlXdq6rHAZQDmNlum5kA/kdb/AnAWSJi35+f3ZCN4xcTnQzbOhq64XXDX14KdQaAd2Oe10SWed1GAbwoIlUiws9RAobjF8ezbRKC2sW1pkPoNmzraOiG1w1/eSnUiSYebd83322b8aqai5aPx/9dRC5PeBCRIhGpFJHK+vp6D2HR6WDj+MUsTm2qws4dZGzMk41tduOUD143/OWlUNcAOC/meSaA9jdrHLdR1dZ/DwJYh5aP0jtQ1TJVDalqaODAgd6iJzoJbsXJNjPKZziuszFPNrbZDfPRxmQuvBTqbQCGi8hQETkDwBwA7W9QbABwU6T39zgAR1S1TkT6ikh/ABCRvgCuBrDzNMZP1GluxYmTMrRxy1NP5dZmGzsa2ngOODGZi17JNlDVRhFZBOAFAKkAHlHVXSJya2T9KgCbAHwRwB4AHwGYF9l9EIB1ItJ6rF+o6vOnvRV00mwcv9gNJ2UgSo7XDX8lLdQAoKqb0FKMY5etivlZAfx7gv32Asg+xRipC9k4frGb9OXpVn0Np7Sg1HQI3UZodYhDZ0bwuuEvjkxmORvHL3YrTrZNQlCU5/xFDBuLuI1tduOUD143/MVCbTkbxy92K062kZJEX9hoYWOebGyzG6d88LrhLxZqso5bcbJxIgonbnnqqdzabGNHQxvPAScmc8FCbTkbxy92w0kZyAk7GrbhdcNfLNSWs6njlBe2TUJQMKLAdAjdRvrydNMhBAavG/5iobacjeMXuxUn2yYhqJhb4bjOxiLu1mbbOhoCzvngdcNfLNSWs3H8YrfiZJvCtYWO62zMk41tduOUD143/MVCTdZxK0622fjWRsd1NubJrc02djS08RxwYjIX1hXqsqqWL+pLiUQfrf8BhWsL45a3bh+7rGJ3BcJHw3HLWu9r5pXlRZe13s8q3lwct21VuApV4aq4Za0fI6UvT48ua/2eYlFFUdy24aNhVOyuiFuWrE1DfjjEl9x2F27FiZMytHHLU0/l1mYbOxraeA44MZoLVQ3cIy8vT7sKitFlrx1U7ducuro2+nNlbaXf4RgXm4/YXKiqbnhzg9/hGOV2btj+u9L+3FiwYYHf4RjnlA/brxuqHc+P03EIp4enIUS7u2HlB7C/oSn6vNca+3os2tjmkzGjfIZVw0Ta1NZTtbp6NYfOJCOsKNT7G5rQOL/lo2gpQfRne8RfjGOLto3jF9vWXjdlVWWOIy7ZmCcb2+zGKR+8bvjLunvUG+a0n6Gz52u9h00tmI82CzcudFxnY55sbLMb5qONyVxYV6jz0u0bTN7tYmwjt3xwUoY2Np43bm22saOhjeeAE5O5sK5QZ6zIMB1CoNg4frEbTspATqrC9vX6dsLrhr+sK9QUj+MXx7NtEgIbbwWdrBnlM0yHEBi8bviLhdoCbhdjG8cvZnFq43YryMY82dhmN0754HXDX9YV6gW5C0yH4Du3i7GN4xfb2E/BidutIBvzZGOb3Tjlg9cNf1lXqG38HiTvy8dzy4eNE1E4sfG8cWuzjR0NbTwHnJjMhXWFunVoTmph4/jFbjgpAzlhR8M2vG74y7pCXV1XbTqEQLFx/GI3tk1CYOOtoJNlW0dDN7xu+Mu6Qm0jt4tx64QiNnHLh22TELjdCrKxiNvYZjdO+eB1w1/WFeq0fmmmQ/Cd28V4dfVqHyMJBhv7KThxuxVkY55sbLMbp3zwuuEv6wp1+E77Jqfgffl4zEcbt1tBNubJrc02djS08RxwYjIX1hXq1rmfbcL78vHc8mHbRANubDxv3NpsY0dDG88BJyZz4alQi8hUEdktIntE5J4E60VEHo6s3yEiuV739VvJlhLTIQSKjeMXu7FtEgIbbwWdLNs6GrrhdcNfSQu1iKQCWAlgGoBRAOaKyKh2m00DMDzyKALw007sa4f6emDbtpZ/feZ2MTYyfrHBXADu+TAy8L7BfLjdCjJWxAP6u2Kko2FAf1eMjXse0HOjq3l5R50PYI+q7lXV4wDKAcxst81MAP+jLf4E4CwRSfO4b8+3di1wwQXAVVe1/Lt2ra+Hd7sY+z5+seFcAAHrp2A4H263gozkKcC/K74L8O+KkXHPLT43RNX9npyIzAIwVVXnR55/FcDFqrooZpuNAB5U1a2R578B8C0AQ5Ltm0j//v01L+/03bjfUnccE9LOAAAcPX4U/c/of9peO6kTJ4A//Qlobm5blpICjBsH9O7tSwj7Du/DkLOGRJ/H5mPLvi2YMGSCL3EEIRdAfD5icwHYl4/27Y3NR/vzpssFIB88N+I55cPXXACByIfbdfR02Lx5s/MX9VXV9QHgOgBrYp5/FcCP223zKwCXxjz/DYA8L/vGrCsCUBl57EwWV7d5ABcpcFgBLQVUWx6HFbjIeGzMBfMRpAfzwVwwHwkfXt5RXwKgWFWnRJ7fGynwP4jZphTAZlVdG3m+G8BEtLyjdt3XJiJSqaoh03EEAXMRj/mIx3y0YS7i2ZgPL/eotwEYLiJDReQMAHMAtJ/vawOAmyK9v8cBOKKqdR73JSIiIge9km2gqo0isgjACwBSATyiqrtE5NbI+lUANgH4IoA9AD4CMM9t3y5pCRERUQ+UtFADgKpuQksxjl22KuZnBfDvXve1mF1f0nXHXMRjPuIxH22Yi3jW5SPpPWoiIiIyx7ohRImIiLoTFmofiMgjInJQRHaajsU0ETlPRF4WkTdEZJeI/F/TMZkkIn1E5M8i8lokH9aPcSsiqSLyl8j4DFYTkX0i8lcR2S4ilabjMUlEzhKRZ0Tkzcj14xLTMfmFH337QEQuB9CAltHbxpiOx6TIiHVpqlotIv0BVAG4RlVfNxyaESIiAPqqaoOI9AawFcD/1ZYR/qwkIosBhAB8WlXtm7IqhojsAxBS1X+YjsU0EfkZgN+p6prIt4g+paqHTcflB76j9oGqvgLgfdNxBIGq1qlqdeTnowDeAJBhNipztEVD5GnvyMPav55FJBPAdABrTMdCwSEinwZwOYD/BgBVPW5LkQZYqMkgERkC4EIAr5qNxKzIR73bARwE8GtVtTkfPwTwTQDNyTa0hAJ4UUSqRKTIdDAGfRZAPYBHI7dF1ohIX9NB+YWFmowQkX4AngXwdVX90HQ8Jqlqk6rmAMgEkC8iVt4eEZECAAdV1dDUTIE0XlVz0TID4b9HbqPZqBeAXAA/VdULAfwTgPFpk/3CQk2+i9yLfRbAE6r6S9PxBEXko7zNAKYaDsWU8QBmRO7LlgOYJCKPmw3JLFUNR/49CGAdWmYktFENgJqYT5ueQUvhtgILNfkq0nnqvwG8oaorTMdjmogMFJGzIj9/EsBkAG+ajcoMVb1XVTNVdQhahhv+rareaDgsY0Skb6TDJSIf814NwMpvjqjqewDeFZHPRxZdCcCaDqieRiajUyMia9EySck5IlIDYImq/rfZqIwZj5ZZ1P4auS8LAN+OjGBnozQAPxORVLT84fyUqlr/tSQCAAwCsK7lb1v0AvALVX3ebEhG3Q7giUiP772IDFVtA349i4iIKMD40TcREVGAsVATEREFGAs1ERFRgLFQExERBRgLNRERUYCxUBMREQUYCzUREVGAsVATEREF2P8H0xwL+Hk7H+cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# nbi:hide_in\n",
    "# library\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "plt.rcParams[\"figure.figsize\"] = (8, 5)\n",
    "plt.gca().spines['top'].set_visible(False)\n",
    "plt.gca().spines['right'].set_visible(False)\n",
    "plt.gca().spines['bottom'].set_visible(False)\n",
    "\n",
    "n = 6\n",
    "num_samples=1000\n",
    "\n",
    "range_x = np.arange(1,n+1)\n",
    "pmf_values = np.array([(2*i-1)/n**2 for i in range(1,n+1)])\n",
    "\n",
    "# generate data\n",
    "toss = np.random.randint(1,n+1,(2,num_samples))\n",
    "data = np.amax(toss, axis=0).squeeze()\n",
    "\n",
    "# compute empirical pmf\n",
    "def epmf(data):\n",
    "    erange_x, counts = np.unique(data, return_counts=True)\n",
    "    epmf_values = counts/data.size\n",
    "    return epmf_values, erange_x\n",
    "\n",
    "epmf_values, erange_x = epmf(data)\n",
    "\n",
    "# plot \n",
    "plt.ylim(0,0.4) \n",
    "plt.axhline(y=0, color='k')\n",
    "plt.xticks(range_x)\n",
    "\n",
    "plt.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20)\n",
    "plt.bar(range_x, pmf_values, width=1, color=(1, 1, 1, 0), \n",
    "        edgecolor='#039be5', linewidth=1, label=\"True histogran\")\n",
    "plt.bar(erange_x, epmf_values, width=0.9, color=(1, 1, 1, 0), \n",
    "        edgecolor='green', linewidth=1,linestyle=\"--\", label=\"Relative frequency histogram\")\n",
    "plt.legend()\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Empirical cdf\n",
    "\n",
    "Similar to empirical pmf, the empirical cfd is computed as\n",
    "\n",
    "$$F_{\\text{emp}}(x)=\\frac{\\text{number of outomes }\\leq x}{\\text{size of data}}.$$\n",
    "\n",
    "We consider the double die toss experiment. Let the random variable $X$ be the maximum of two tosses. Again,   $\\text{Range}(X)=\\{1,\\dots, 6\\}$ and it can be checked that\n",
    "     $$F(x)=\n",
    "\\begin{cases}\n",
    "0 & x<1\\\\\n",
    "\\sum\\limits_{i=1}^k\\frac{2i-1}{6^2}=\\frac{k^2}{6^2} & k\\leq x<k+1, \\text{ for }k=1,\\dots,5\\\\\n",
    "1 & 6<x.\n",
    "\\end{cases}$$\n",
    "In the numerical example below we collect 1000 pairs of tosses of the die and collect the maximum value of each pair as our data. The graph shows the real and the empirical cdf-s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAEvCAYAAAB2Xan3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAfDElEQVR4nO3df3DU9b3v8dc7AUQRyw+pJsEiOlql/EiTFXEqhauUgiegduwIR2vtXFjsaGvFjtrb3mY5TudWT3Hac7VAoKIde8HWWgsU7S1VbOspV7OIInigVJwaEjUiVPHHYMj7/pF1WUMgy7rsZ7+7z8fMTvP9fD/ZfSXYeeX73e9+vubuAgAA4VSEDgAAQLmjjAEACIwyBgAgMMoYAIDAKGMAAAKjjAEACKxPqBeeNm2aP/bYY6FeHgCAEKynwWBHxm+88UaolwYAoKhwmhoAgMAoYwAAAqOMAQAILNgFXD354IMP1NLSovfffz90FEjq37+/hg8frr59+4aOAgAlrajKuKWlRQMHDtTpp58usx4vOEOBuLt2796tlpYWjRw5MnQcAChpRXWa+v3339fQoUMp4iJgZho6dChnKQCgAIqqjCVRxEWEfwsAKIyiK+OQdu/erdraWtXW1urUU09VTU1Nenv//v0FzbJu3TpddtllkqT33ntPF110kWpra/XQQw8VNAcA4NgrqveMQxs6dKg2bdokSUokEjrxxBP17W9/+yNz3F3uroqKwv0dk0wmZWbpbACA0sKRcRZ27Nih0aNH67rrrlNdXZ1eeeUVDRo0KL1/5cqVmjNnjiTptdde05e+9CXFYjGNHz9eGzZsOOT5Ojo6dNNNN2n06NEaO3asfvrTn0qSfve73+nTn/60LrzwQv32t7+VJLW2turaa69Vc3Ozamtr9fLLLx/7HxgAUFAcGWdp69atWr58uRYvXqyOjo7DzvvmN7+pW265RRMmTNDLL7+shoYGvfDCCx+Zs2jRIrW2tuq5555TZWWl3nzzTb377ruaN2+ennzySZ1xxhm64oorJEnV1dVavHix7r77bj3yyCPH9GcEAIRBGWfpzDPP1HnnndfrvHXr1mnbtm3p7T179ui9997T8ccf/5E53/rWt1RZWSlJGjJkiJqbm3X22WfrzDPPlCRdddVV+vnPf57nnwIAUIyKuoz7LGvN+3N2zKnO6fsGDBiQ/rqiokLunt7O/PiPu+vpp59Wv379Dvtc7t7jlcpcvQwA5amoyzjX4jzWKioqNHjwYP3tb3/TmWeeqd/85jcaNmyYJGnKlCm65557dNNNN0mSNm3apNra2o98/9SpU7Vo0SJNnDgxfZp61KhR2r59u3bu3KnTTz9dK1asKPjPBQAIgwu4cnTHHXdo2rRpuvjiizV8+PD0+D333KOnnnpKY8eO1ahRo7R06dJDvnfevHk69dRTNXbsWI0bN06//OUvdcIJJ2jx4sWaPn26Jk6cqDPOOKOQPw4AICDLPN1aSLFYzJubmz8y9uKLL+rcc88Nkgc9498EAPKqx/cjOTIGACAwyhgAgMAoYwAAAuu1jM3sXjN73cxeOMx+M7P/MLMdZva8mdXlPyYAAKUrmyPj+yRNO8L+6ZLOSj3ikhZ9/FgAAJSPXsvY3f8k6c0jTLlU0s+9ywZJg8ysKl8BAQAodfl4z7hG0isZ2y2psUiqrKxM3zaxtrZWP/zhD/PyvJdccon27t171PuO5L777tMNN9yQU57M721vb9f555+vz372s/rzn/+c0/MBAHKXjxW4evrMVI8fXjazuLpOZetTn/pUHl46/44//vhjcqvCtWvXHjL24e0Ye9pXSH/84x91zjnn6P777w+aA8DHU3H7YHnnXqniFFVWbVTnWwvlb991cP+wRyVJne3T02M2cL4qTrpZB9rqpM7Xugb7jlHlJx9T555b5O/+4uD3n5qU9j+vzje/dvD7B92higFX68CujGOw/lNUOfR+Hdj9Ven9denhyppd6nznAfneWw8+55DlUr+x6ny1/uBznnCVKgbfqQOvT5M+2JyaWNw/kw68ps75C5WzDwvhSA9Jp0t64TD7lkianbG9TVJVb89ZX1/v3W3duvWQsUIbMGBAj+MjRozw73znOz5hwgSvr6/3ZDLpU6dO9TPOOMMXLVrk7u5PPPGET5w40S+77DI/99xzfd68eX7gwIH097e3t/vOnTv9nHPO8a9//eteW1vrL7/8cnqfu/v999/vY8aM8bFjx/rVV1/t7u6rVq3y8ePHe21trV988cX+6quvurv78uXL/frrrz8k69tvv+3XXnutjx492seMGeMPPfSQu7vfe++9ftZZZ/nnP/95nzNnjl9//fX+7LPP+mmnneYnn3yyjxs3zt99992PPFcx/JsAyI5+MCZ0BPSux07Mx2nqVZKuSV1VPUHSP929LQ/PG8R77733kdPUDz74YHrfaaedpr/+9a+aOHGirr32Wj300EPasGGDvv/976fnPP3001q4cKE2b96sv//973r44YcPeY1t27bpmmuu0bPPPqsRI0akx7ds2aIf/OAHevzxx/Xcc8/pJz/5iSTpwgsv1IYNG/Tss89q1qxZuvPOO4/4M9x+++36xCc+oc2bN+v555/XRRddpLa2NjU2Nuqpp57SH/7wB23dulWSVFtbq3/7t3/TlVdeqU2bNn3k7lIAoqXyk4+FjoAc9Xqa2sxWSJos6WQza5HUKKmvJLn7YklrJV0iaYekdyV9rednioYjnaaeOXOmJGnMmDHat2+fBg4cqIEDB6p///7p93zHjx+fXld69uzZ+stf/pK+N/GHRowYoQkTJhzy/I8//riuuOIKnXzyyZK6bq0oSS0tLbryyivV1tam/fv3a+TIkUf8GdatW6eVK1emtwcPHqxHHnlEkydPTt/Q4sorr9T27dt7/X0AiI7OPbdIeiB0DOSg1zJ299m97HdJ1+ctUYbE+oQWPLkgvd08t2st69jSWHqscVKjEpMTql5YrbZ9XQfkdVV1SsaTiq+Oa+nGgzdq2DV/l6oH5n4nqOOOO05S112bPvz6w+2Ojg5Jh94GsafbImbejjGTH+bWit/4xjc0f/58zZw5U+vXr1cikThizsM9D7doBEpb13uhlHEUFfUtFBOTE0pMThwy7o2HXh/WevOh9z5umtGkphlNxyLaYT399NPauXOnRowYoQcffFDxeDzr77344ot1+eWX66abbtLQoUP15ptvasiQIfrnP/+pmpquCwmyuchq6tSpuvvuu/XjH/9YkrRnzx6df/75uvHGG7V7926ddNJJ+tWvfqVx48bl9kMCAPKK5TC76f6e8W233XZU33/BBRfotttu0+jRozVy5EhdfvnlWX/vZz7zGX33u9/VpEmTNG7cOM2fP1+SlEgk9OUvf1kTJ05Mn8I+ku9973vas2ePRo8erXHjxumJJ55QVVWVEomELrjgAk2ZMkV1dSyUBgDFglso5tH69ev1ox/9SGvWrAkdJW+i/m8ClJPKJRt1YB5/aBc5bqEIACVt//OhEyBHlHEeTZ48uaSOigFES+bCFYgWyhgAgMCKroxDvYeNQ/FvAQCFUVRl3L9/f+3evZsSKALurt27d6t///6howDIkg26I3QE5KioPmc8fPhwtbS0qL29PXQUqOuPo+HDh4eOASBLFQOuDh0BOSqqMu7bt2+vSz0CAHrWdZchzixGUVGdpgYAoBxRxgAABEYZA0Cp6D8ldALkiDIGgBJRObT3G8mgOFHGAFAiDuz+augIyBFlDACl4v11oRMgR5QxAACBUcYAAARGGQNAiais2RU6AnJEGQNAieh854HQEZAjyhgASoTvvTV0BOSIMgYAIDDKGACAwChjACgRFUOWh46AHFHGAFAq+o0NnQA5oowBoER0vlofOgJyRBkDABAYZQwAQGCUMQCUCDvhqtARkCPKGABKRMXgO0NHQI4oYwAoEQdenxY6AnJEGQNAqfhgc+gEyBFlDABAYJQxAJSKilNCJ0COKGMAKBGVVRtDR0COKGMAKBGdby0MHQE5oowBoET423eFjoAcUcYAAARGGQMAEBhlDAAlomLYo6EjIEdZlbGZTTOzbWa2w8xu62H/J8xstZk9Z2ZbzOxr+Y8KAEBp6rWMzaxS0j2SpksaJWm2mY3qNu16SVvdfZykyZIWmlm/PGcFABxBZ/v00BGQo2yOjMdL2uHuL7n7fkkrJV3abY5LGmhmJulESW9K6shrUgAASlQ2ZVwj6ZWM7ZbUWKa7JZ0rqVXSZkk3untnXhICAFDisilj62HMu21/UdImSdWSaiXdbWYnHfJEZnEzazaz5vb29qMOCwA4PBs4P3QE5CibMm6RdFrG9nB1HQFn+pqkh73LDkk7JZ3T/YncvcndY+4eGzZsWK6ZAQA9qDjp5tARkKNsyvgZSWeZ2cjURVmzJK3qNucfki6WJDM7RdKnJb2Uz6AAgCM70FYXOgJy1Ke3Ce7eYWY3SPq9pEpJ97r7FjO7LrV/saTbJd1nZpvVdVr7Vnd/4xjmBgB01/la6ATIUa9lLEnuvlbS2m5jizO+bpU0Nb/RAAAoD1mVMQCUiyF3DNGe9/dIkprnNkuSYktj6f2NkxqVmJxQ9cJqte1rkyTVVdUpGU8qvjqupRuXpufumr9LydakZq6cmR5b0rBE8fq4bMHBa2Mbzm7Q6tmrNWPFDK3ZviY97o2upmST5q2Zlx5bNWuV6qvrVXPXwQ+1zK2bq6YZTVLl6Xn6LaDQzL37hdGFEYvFvLm5OchrA8DhJNYnlJicCB0jJ32WtapjTnXoGDiynj6hxNrUAJApqkWMaKOMASBD9UKOLFF4lDEAZPjwfWCgkChjAAACo4wBIENdFQtnoPAoYwDIkIwnQ0dAGaKMASBDfHU8dASUIcoYADJkLtoBFAplDABAYJQxAACBUcYAkGHX/F2hI6AMUcYAkCHZytXUKDzKGAAyZN5hCSgUyhgAgMAoYwAAAqOMASDDkoYloSOgDFHGAJAhXs8KXCg8yhgAMtgCCx0BZYgyBgAgsD6hAwBAsemzrDV0hJwMPo6j+qiijAEgU/8p6phTHToFygynqQEgQ+XQ+0NHQBmijAEgw4HdXw0dAWWIMgaATO+vC50AZYgyBgAgMMoYAIDAKGMAyFBZw/2MUXiUMQBk6HzngdARUIYoYwDI4HtvDR0BZYgyBgAgMMoYAIDAKGMAyFAxZHnoCChDlDEAZOo3NnQClCHKGAAydL5aHzoCyhBlDABAYJQxAACBUcYAkMFOuCp0BJQhyhgAMlQMvjN0BJQhyhgAMhx4fVroCChDWZWxmU0zs21mtsPMbjvMnMlmtsnMtpjZk/mNCQAF8sHm0AlQhvr0NsHMKiXdI+kLklokPWNmq9x9a8acQZJ+Kmmau//DzD55rAIDAFBqsjkyHi9ph7u/5O77Ja2UdGm3Of8q6WF3/4ckufvr+Y0JAAVScUroBChD2ZRxjaRXMrZbUmOZzpY02MzWm1nSzK7JV0AAKKTKqo2hI6AMZVPG1sOYd9vuI6le0r9I+qKk/2lmZx/yRGZxM2s2s+b29vajDgsAx1rnWwtDR0AZyqaMWySdlrE9XFJrD3Mec/d33P0NSX+SNK77E7l7k7vH3D02bNiwXDMDwDHjb98VOgLKUDZl/Iyks8xspJn1kzRL0qpuc34raaKZ9TGzEySdL+nF/EYFAKA09Xo1tbt3mNkNkn4vqVLSve6+xcyuS+1f7O4vmtljkp6X1Clpmbu/cCyDAwBQKnotY0ly97WS1nYbW9xt+98l/Xv+ogGIquqF1Wrb1yZJqquqUzKeVHx1XEs3Lk3P2TV/l5KtSc1cOTM9tqRhieL1cdmCg5eqNJzdoNWzV2vGihlas31NetwbXU3JJs1bMy89tmrWKtVX16vmroPXmM6tm6umGU2qb6rXxraui7OqTqxS682tSqxPaMGTC9Jzm+c2q2LYo3n8TQDZMffu12IVRiwW8+bm5iCvDeDYSaxPKDE5ETpGzvosa1XHnOrQMVC6eroomuUwAeRX5pEmgOxQxgAABEYZAwAQGGUMIK+a53ItCHC0KGMAAAKjjAHkVWxpLHQEIHIoYwAAAqOMAQAIjDIGkFeNkxpDRwAihzIGkFdRXn0LCIUyBpBX1QtZShI4WpQxgLz68AYRALJHGQMAEBhlDCCv6qrqQkcAIocyBpBXyXgydAQgcihjAHkVXx0PHQGIHMoYQF4t3bg0dAQgcihjAAACo4wBAAiMMgaQV7vm7wodAYgcyhhAXiVbuZoaOFqUMYC8mrlyZugIQORQxgAABEYZAwAQGGUMIK+WNCwJHQGIHMoYQF7F61mBCzhalDGAvLIFFjoCEDmUMQAAgVHGAAAERhkDyKuGsxtCRwAihzIGkFerZ68OHQGIHMoYQF7NWDEjdAQgcvqEDgCgtKzZvkZ9lrWGjpGzwcdxNTgKjzIGkHcdc6pDRwAihdPUAAAERhkDyKvKGu5nDBwtyhhAXnW+80DoCEDkUMYA8sr33ho6AhA5lDEAAIFlVcZmNs3MtpnZDjO77QjzzjOzA2Z2Rf4iAgBQ2notYzOrlHSPpOmSRkmabWajDjPvDkm/z3dIANFRMWR56AhA5GRzZDxe0g53f8nd90taKenSHuZ9Q9KvJb2ex3wAoqbf2NAJgMjJpoxrJL2Ssd2SGkszsxpJl0tanL9oAKKo89X60BGAyMmmjHtaG867bf9Y0q3ufuCIT2QWN7NmM2tub2/PNiMAACUtm+UwWySdlrE9XFL3hWdjklaamSSdLOkSM+tw90cyJ7l7k6QmSYrFYt0LHQCAspRNGT8j6SwzGylpl6RZkv41c4K7j/zwazO7T9Ka7kUMoDzYCVeFjgBETq9l7O4dZnaDuq6SrpR0r7tvMbPrUvt5nxhAWsXgO0NHACInq7s2uftaSWu7jfVYwu5+7cePBSCqDrw+TdLzoWMAkcIKXADy64PNoRMAkUMZAwAQGGUMIL8qTgmdAIgcyhhAXlVWbQwdAYgcyhhAXnW+tTB0BCByKGMAeeVv3xU6AhA5lDEAAIFRxgAABEYZA8irimGPho4ARA5lDABAYJQxgLzqbJ8eOgIQOZQxAACBUcYAAARGGQPIKxs4P3QEIHIoYwB5VXHSzaEjAJFDGQPIqwNtdaEjAJFDGQPIr87XQicAIocyBgAgMMoYQH71HRM6ARA5lDGAvKr85GOhIwCRQxkDyKvOPbeEjgBEDmUMIK/83V+EjgBETp/QAQAcasgdQ9S/T3+13tyqxPqEFjy5IL2veW6zJCm2NJYea5zUqMTkhKoXVqttX5skqa6qTsl4UvHVcS3duDQ9d9f8XUq2JjVz5cz02JKGJYrXx2ULLD3WcHaDVs9erRkrZmjN9jXpcW90NSWbNG/NvPTYqlmrVF9dr5q7aiQblMffBFAezN2DvHAsFvPm5uYgrw0UO1tg8sYw/9/8uPosa1XHnOrQMYBiZT0NcpoaKEJ1VSycAZQTyhgoQsl4MnQEAAVEGQNFKL46HjoCgAKijIEilHnBFYDSRxkDABAYZQwAQGCUMVCEds3fFToCgAKijIEilGzlamqgnLACF1CEZq6cqcqaaB4dDz6uxzUNABwBZQwUKVaxAsoHp6kBAAiMMgaKkA26I3QEAAVEGQNFqGLA1aEjACggyhgoQgd21YSOAKCAKGMAAAKjjAEACIwyBopR/ymhEwAooKzK2Mymmdk2M9thZrf1sP8qM3s+9fhPMxuX/6hA+agcen/oCAAKqNcyNrNKSfdImi5plKTZZjaq27Sdkia5+1hJt0tqyndQoJwc2P3V0BEAFFA2R8bjJe1w95fcfb+klZIuzZzg7v/p7ntSmxskDc9vTKDMvL8udAIABZRNGddIeiVjuyU1djj/XdKjPe0ws7iZNZtZc3t7e/YpAQAoYdmUcU+rvnuPE83+m7rK+Nae9rt7k7vH3D02bNiw7FMCAFDCsrlRRIuk0zK2h0tq7T7JzMZKWiZpurvvzk88oDxF9Y5NAHKTzZHxM5LOMrORZtZP0ixJqzInmNmnJD0s6Svuvj3/MYHy0vnOA6EjACigXsvY3Tsk3SDp95JelPRLd99iZteZ2XWpad+XNFTST81sk5k1H7PEQBnwvT2+0wOgRGV1P2N3XytpbbexxRlfz5E0J7/RAAAoD6zABQBAYJQxUIQqhiwPHQFAAVHGQDHqNzZ0AgAFRBkDRajz1frQEQAUEGUMAEBglDEAAIFRxkARshOuCh0BQAFl9TljIGpa325VsjWpmStnpseWNCxRvD4uW3BwufWGsxu0evZqzVgxQ2u2r0mPe6OrKdmkeWvmpcdWzVql+up61dx18D4pc+vmqmlGk+qb6rWxbaMkqerEKrXe3KrE+oQWPLkgPbd5btdaOLGlsfRY46RGJSYnVL2wWm372iRJdVV1qhi8Ol+/CgARYO493vPhmIvFYt7czEJdODZsgckbw/y3nQ99lrWqY0516BgA8q+nmy9xmhoAgNAoYwAAAqOMUZLm1s0NHQEAskYZoyQ1zWgKHQEAskYZoyTVN7GCFYDooIxRkj78mBEARAFlDABAYCz6gZJUUXmK+ixrDR0jZ4OP6/GjiABKFGWMkmSnbmTRDACRwWlqlKTOtxaGjgAAWaOMUZL87btCRwCArFHGAAAERhkDABAYZYySVDHs0dARACBrlDEAAIFRxihJne3TQ0cAgKxRxgAABEYZAwAQGGWMkmQD54eOAABZo4xRkipOujl0BADIGmtT47CSrUlJUmxpLD3WOKlRickJVS+sVtu+NklSXVWdkvGk4qvjWrpxaXrurvm7lGxNaubKmemxJQ1LFK+PyxYcvBFCw9kNWj17tWasmKE129ekx73R1ZRs0rw189Jjq2atUn11vWruqkmPza2bq6YZTapvqj9460Q7QdI7+flFAMAxZu4e5IVjsZg3NzcHeW1kxxaYvDHMfx8fV59lrdwoAkAx6vGWbJymBgAgMMoYAIDAKGMcVuOkxtARAKAsUMY4rMTkROgIAFAWKGMcVvVCLoACgEKgjHFYH350CQBwbPE5YxxRn2WtoSPkZPBxPX56AACKEmWMw+s7hs/qAkABcJoah1X5ycdCRwCAspBVGZvZNDPbZmY7zOy2Hvabmf1Hav/zZlaX/6gotM49t4SOAABlodcyNrNKSfdImi5plKTZZjaq27Tpks5KPeKSFuU5JwqtvV3+7i+k9vbQSY5ee7v0zDPRzC5FO3+Us0vRzh/l7MjqyHi8pB3u/pK775e0UtKl3eZcKunn3mWDpEFmVpXnrCiUFSukESO6vh4xoms7Kj7M/oUvRC+7FO38Uc4uRTt/lLNDUhY3ijCzKyRNc/c5qe2vSDrf3W/ImLNG0g/d/S+p7T9KutXdD3sniIEDB3p9fX0efgTk1QcfSBs2SJ2dejIhTUpIqqiQJkyQ+vYNHK4XGdnTopJdinb+KGeXop0/ytnL0Pr163v+qIe7H/Eh6cuSlmVsf0XS/+4253eSLszY/qOk+h6eKy6pOfV4obfX5hHgIZ3n0l6XfInk3vXY69J5wbOVcvao549y9qjnj3J2HulHNkfGF0hKuPsXU9vfSZX4/8qYs0TSendfkdreJmmyu7NqRISZWbO7x3qfWXyinF2Kdv4oZ5einT/K2ctdNu8ZPyPpLDMbaWb9JM2StKrbnFWSrkldVT1B0j8pYgAAstProh/u3mFmN0j6vaRKSfe6+xYzuy61f7GktZIukbRD0ruSvnbsIgMAUFqyWoHL3deqq3AzxxZnfO2Srs9vNBSBptABPoYoZ5einT/K2aVo549y9rLW63vGAADg2GI5TAAAAqOMcQgzu9fMXjezF0JnOVpmdpqZPWFmL5rZFjO7MXSmbJlZfzN72syeS2VfEDrT0TKzSjN7NrX2QKSY2ctmttnMNpnZYddIKFZmNsjMHjKz/0r9939B6EzIHqepcQgz+7ykfepaVW106DxHI7XyW5W7bzSzgZKSki5z962Bo/XKzEzSAHffZ2Z9Jf1F0o3etapdJJjZfEkxSSe5e0PoPEfDzF6WFHP3N0JnyYWZ3S/pz+6+LPXJlxPcfW/oXMgOR8Y4hLv/SdKboXPkwt3b3H1j6uu3Jb0oqSZsqux4l32pzb6pR2T+Wjaz4ZL+RdKy0FnKjZmdJOnzkn4mSe6+nyKOFsoYJcvMTpf0WUn/L2yS7KVO826S9LqkP7h7ZLJL+rGkWyR19jaxSLmk/2tmSTOLhw5zlM6Q1C5peeptgmVmNiB0KGSPMkZJMrMTJf1a0rfc/a3QebLl7gfcvVbScEnjzSwSbxOYWYOk1909GTrLx/A5d69T113ork+9XRMVfSTVSVrk7p+V9I6kQ253i+JFGaPkpN5v/bWkX7j7w6Hz5CJ1inG9pGmBo2Trc5Jmpt53XSnpIjN7IGyko+Puran/fV3Sb9R1x7qoaJHUknEm5SF1lTMigjJGSUldBPUzSS+6+12h8xwNMxtmZoNSXx8vaYqk/wqbKjvu/h13H+7up6trydzH3f3qwLGyZmYDUhf8KXV6d6qkyHyawN1flfSKmX06NXSxpKK/aBEHZbUCF8qLma2QNFnSyWbWIqnR3X8WNlXWPqeuO4ttTr33Kkn/I7WKXLGrknS/mVWq6w/lX7p75D4iFFGnSPpN199y6iPp/7j7Y2EjHbVvSPpF6krql8SyxJHCR5sAAAiM09QAAARGGQMAEBhlDABAYJQxAACBUcYAAARGGQMAEBhlDABAYJQxAACB/X/PXWah+qIsMAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# nbi:hide_in\n",
    "# library\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "plt.rcParams[\"figure.figsize\"] = (8, 5)\n",
    "plt.gca().spines['top'].set_visible(False)\n",
    "plt.gca().spines['right'].set_visible(False)\n",
    "plt.gca().spines['bottom'].set_visible(False)\n",
    "\n",
    "n = 6\n",
    "num_samples=100\n",
    "range_x = np.arange(1,n+1)\n",
    "\n",
    "# generate data\n",
    "toss = np.random.randint(1,n+1,(2,num_samples))\n",
    "data = np.amax(toss, axis=0).squeeze()\n",
    "\n",
    "# compute true cdf values\n",
    "cdf_values = np.array([i**2/n**2 for i in range(1,n+1)])\n",
    "\n",
    "# compute empirical cdf values\n",
    "def ecdf(data):\n",
    "    erange_x, counts = np.unique(data, return_counts=True)\n",
    "    cdf_emp = np.cumsum(counts)/data.size\n",
    "    return cdf_emp, erange_x\n",
    "\n",
    "ecdf_values, erange_x = ecdf(data)\n",
    "\n",
    "# add padding \n",
    "def padding(cdf_values, range_x):\n",
    "    edge = (range_x[0]-2, range_x[-1]+2)\n",
    "    range_x_padded = np.pad(range_x, (1,1), 'constant', constant_values=edge)\n",
    "    cdf_values_padded = np.pad(cdf_values, (1,1), 'constant', constant_values=(0, 1))\n",
    "    return cdf_values_padded, range_x_padded\n",
    "\n",
    "cdf_values_padded, range_x_padded = padding(cdf_values, range_x)\n",
    "ecdf_values_padded, erange_x_padded = padding(ecdf_values, erange_x)\n",
    "\n",
    "# plot setup\n",
    "plt.ylim(0,1.1)\n",
    "plt.xlim(-1.2,8.2)\n",
    "plt.axhline(y=0, color='k')\n",
    "plt.xticks(range_x)\n",
    "\n",
    "# plot cdf using step function\n",
    "plt.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20)\n",
    "plt.step(range_x_padded, cdf_values_padded, where='post', color='#039be5', linewidth=1, label=\"True cdf\")\n",
    "plt.step(erange_x_padded, ecdf_values_padded, where='post', color=\"green\", linewidth=1, linestyle=\"--\", label='Empirical cdf')\n",
    "plt.legend()\n",
    "\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bernoulli trials\n",
    "\n",
    "\n",
    "**Definition**\n",
    "    \n",
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "In Bernoulli experiment the outcome is either of Type I (success) with probability $p$ or Type II (fail) with probability $q=1-p$. \n",
    "\n",
    "</div>\n",
    "\n",
    "A **Bernoulli trials** are the outcomes of several independent Bernoulli experiments with the same success probabilities.\n",
    "\n",
    "**Example**\n",
    "Consider the coin flip experiment with a biased coin that lands head with probability $p$ and tail with probability $q=1-p$.A single coin flip will be a Bernoulli experiment.$n$ independent coin flips will be a Bernoulli trial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geometric distribution\n",
    "\n",
    "A random experiment has a probability of success equal to $p$ and probability of failure $q=1-p$. Consider the following new experiment: we are doing consecutive  random trials until we reach a success.The set of outcomes has the form $s=\\overbrace{FF\\cdots F}^{n-1} S$ where number of $F$'s can be any number $n=1,2,\\dots$.  Due to independence, $P(s)=q^np$.  \n",
    "\n",
    "Let $X(s)$ denote the number of trials it took to reach success\n",
    "\n",
    "$$X(\\overbrace{FF\\cdots F}^{n-1} S)=n.$$\n",
    "\n",
    "Observe that \n",
    "\n",
    "$$f(n)=q^{n-1}p,\\quad n=1,2,\\dots.$$\n",
    "\n",
    "**Definition (Geometric distribution)**\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "The pmf of the random variable $X$ is called Geometric distribution.\n",
    "\n",
    "</div>    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3c71640c2b24435186e2fef03223e0c3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.4, description='p', max=1.0, min=0.01, readout_format='', step=0.01)…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# nbi:hide_in\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from ipywidgets import interact, FloatSlider\n",
    "\n",
    "\n",
    "def geometric_pmf(p):\n",
    "    q = 1-p\n",
    "    N=15\n",
    "    range_x = np.arange(1, N, 1)\n",
    "\n",
    "    def geo_pmf(n):\n",
    "        pmf_val = q**(n-1)*p\n",
    "        return pmf_val\n",
    "    mean = 1/p\n",
    "\n",
    "    pmf_values = np.array([geo_pmf(n) for n in range_x])\n",
    "\n",
    "    # plot setup\n",
    "    plt.figure(figsize=(10,7)) \n",
    "    plt.axhline(y=0, color='k')\n",
    "    plt.ylim(-0.01, max(np.max(pmf_values)+0.05, 0.2))\n",
    "    plt.xlim(0, N+1)\n",
    "    plt.xticks(np.arange(0, N+1, 5))\n",
    "    plt.gca().spines['top'].set_visible(False)\n",
    "    plt.gca().spines['right'].set_visible(False)\n",
    "    plt.gca().spines['bottom'].set_visible(False)\n",
    "\n",
    "    # PLotting with plt.bar instead of plt.hist works better when f(x) are knowwn\n",
    "    plt.scatter(np.array([mean]),np.zeros(1), color =\"red\", s=200, label=\"Mean\", zorder=3)\n",
    "    plt.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20)\n",
    "    plt.bar(range_x, pmf_values, width=1, color='#039be5', edgecolor=\"w\", linewidth=1.3)\n",
    "    plt.title(\"Geometric distribution\")\n",
    "    plt.figtext(0.8,0.8, \"p={}\".format(p), ha=\"left\", va=\"top\",\n",
    "            backgroundcolor=(0.1, 0.1, 1, 0.15), fontsize=\"large\")\n",
    "    plt.plot();\n",
    "\n",
    "# create interactive variables\n",
    "p = FloatSlider(min=0.01, max=1, step=0.01, value=0.4, readout_format='')\n",
    "\n",
    "# display the interactive plot\n",
    "interact(geometric_pmf, p=p);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Binomial distribution\n",
    "\n",
    "\n",
    "**Definition**\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "Let $X$ denote the number of successes in $n$ Bernoulli trials with success probability $p$. Then we say that $X$ has the \\textbf{binomial distribution} with parameters $(n,p)$. \n",
    "\n",
    "</div>\n",
    "    \n",
    "If $X$ has binomial distribution with parameters $(n,p)$ then \n",
    "* $\\text{range}(X)=\\{0,1,\\dots, n\\}$.\n",
    "* $$f(x)={n\\choose x}p^xq^{(n-x)}$$ ${n\\choose x}$ ways there can be $x$ successes: the rest $n-x$ are failures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b443244646984ab49f42dbee61b29ab3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=20, description='n', min=1, readout_format=''), FloatSlider(value=0.3, d…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# nbi:hide_in\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from ipywidgets import interact, FloatSlider, IntSlider\n",
    "from scipy.special import comb\n",
    "\n",
    "def binomial_pmf(n, p):\n",
    "    q = 1-p\n",
    "    range_x = np.arange(0, n, 1)\n",
    "\n",
    "    def binom_pmf(x):\n",
    "        pmf_val = comb(n, x, exact=True) * p**x * q**(n-x)\n",
    "        return pmf_val\n",
    "    mean = n*p\n",
    "\n",
    "    pmf_values = np.array([binom_pmf(x) for x in range_x])\n",
    "    \n",
    "    # plot setup\n",
    "    plt.figure(figsize=(10,7)) \n",
    "    plt.axhline(y=0, color='k')\n",
    "    plt.ylim(-0.01, max(np.max(pmf_values)+0.05, 0.2))\n",
    "    plt.xlim(-1, n+1)\n",
    "    plt.xticks(np.arange(0, n+1, 5))\n",
    "    plt.gca().spines['top'].set_visible(False)\n",
    "    plt.gca().spines['right'].set_visible(False)\n",
    "    plt.gca().spines['bottom'].set_visible(False)\n",
    "    plt.title(\"Binomial distribution\")\n",
    "\n",
    "    # PLotting with plt.bar instead of plt.hist works better when f(x) are knowwn\n",
    "    plt.scatter(np.array([mean]),np.zeros(1), color =\"red\", s=200, zorder=3)\n",
    "    plt.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20, zorder=2)\n",
    "    plt.bar(range_x, pmf_values, width=1, color='#039be5', edgecolor=\"w\", linewidth=1.3, zorder=1)\n",
    "    #plt.title(\"Binomial distribution\")\n",
    "    plt.figtext(0.8,0.8, \" n={} \\n p={}\".format(n, p), ha=\"left\", va=\"top\",\n",
    "                backgroundcolor=(0.1, 0.1, 1, 0.15), fontsize=\"large\")\n",
    "    plt.plot();\n",
    "\n",
    "# create interactive variables\n",
    "p = FloatSlider(min=0.01, max=1, step=0.01, value=0.3, readout_format='')\n",
    "n = IntSlider(min=1, max=100, step=1, value=20, readout_format='')\n",
    "# display the interactive plot\n",
    "interact(binomial_pmf, n=n, p=p);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### cdf of the binomial distribution\n",
    "\n",
    "The cdf of the binomial distribution can be computed as \n",
    "$$F(x)=E[X\\leq x]=\\sum_{y=0}^{\\lfloor x\\rfloor}{n\\choose x}p^xq^{(n-x)}\\text{ for } x\\in [0,n].$$\n",
    "\n",
    "However, $F(x)$ does not have a simple closed-form formula. In old times people would use statistical tables that show  precomputed values of the cdf for different choices of parameters $(n,p)$. Nowdays, computers can give you very precise values for the cdf for any $n,p,x$.\n",
    "\n",
    "In the plot below, the shaded area is the value of the cdf at point $x$. The point $x$ is the black dot on the plot. The value of the cdf up to 4 digit precision is given at the top of the plot. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4058198aea7d498696df1025106e2544",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=20, description='n', min=1, readout_format=''), FloatSlider(value=0.3, d…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# nbi:hide_in\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from ipywidgets import interact, FloatSlider, IntSlider\n",
    "from scipy.special import comb\n",
    "from scipy.stats import binom\n",
    "\n",
    "def binomial_cdf(n, p, x):\n",
    "    q = 1-p\n",
    "    range_x = np.arange(0, n, 1)\n",
    "\n",
    "    def binom_pmf(x):\n",
    "        pmf_val = comb(n, x, exact=True) * p**x * q**(n-x)\n",
    "        return pmf_val\n",
    "    mean = n*p\n",
    "\n",
    "    pmf_values = np.array([binom_pmf(x) for x in range_x])\n",
    "    \n",
    "    # plot setup\n",
    "    plt.figure(figsize=(10,7)) \n",
    "    plt.axhline(y=0, color='k')\n",
    "    plt.ylim(-0.01, max(np.max(pmf_values)+0.05, 0.2))\n",
    "    plt.xlim(-4, n+2)\n",
    "    plt.gca().spines['top'].set_visible(False)\n",
    "    plt.gca().spines['right'].set_visible(False)\n",
    "    plt.gca().spines['bottom'].set_visible(False)\n",
    "\n",
    "    # PLotting with plt.bar instead of plt.hist works better when f(x) are knowwn\n",
    "#     plt.scatter(x,np.zeros(1), color =\"red\", s=100, zorder=3)\n",
    "    plt.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20,zorder=2)\n",
    "    barlist = plt.bar(range_x, pmf_values, width=1, color=(1, 1, 1, 0), edgecolor='#039be5',\n",
    "                      linewidth=1, label=\"Histogran\", zorder=1)\n",
    "    for ind, val in enumerate(range_x):\n",
    "        if val<=x:\n",
    "            barlist[ind].set_color('#039be5')\n",
    "            barlist[ind].set_edgecolor(\"w\")    \n",
    "    cdf = binom.cdf(x, n, p)\n",
    "    plt.title(\"cdf and pdf of binomial distribution\" )\n",
    "    plt.figtext(0.7,0.8, \" n={} \\n p={} \\n x={} \\n cdf = {:.4f}\".format(n, p, x, cdf), ha=\"left\", va=\"top\",\n",
    "                backgroundcolor=(0.1, 0.1, 1, 0.15), fontsize=\"large\")\n",
    "    plt.plot();\n",
    "\n",
    "# create interactive variables\n",
    "p = FloatSlider(min=0.01, max=1, step=0.01, value=0.3, readout_format='')\n",
    "n = IntSlider(min=1, max=100, step=1, value=20, readout_format='')\n",
    "x = FloatSlider(min=-4, max=10, step=0.1, value=6.75, readout_format='')\n",
    "\n",
    "# enforce values for x\n",
    "def update_range(*args):\n",
    "    x.max = max(n.value+4,10)\n",
    "n.observe(update_range, 'value')\n",
    "\n",
    "# display the interactive plot\n",
    "interact(binomial_cdf, n=n, p=p, x=x);"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}