{ "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": [ "
" ] }, "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": [ "
" ] }, "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": [ "
" ] }, "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" ] }, "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", "
\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", "
\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", "
\n", "\n", "The pmf of the random variable $X$ is called Geometric distribution.\n", "\n", "
" ] }, { "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", "
\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", "
\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 }