diff --git a/session-2-intuition/bootstrap.ipynb b/session-2-intuition/bootstrap.ipynb index ebd4daf..47748c3 100644 --- a/session-2-intuition/bootstrap.ipynb +++ b/session-2-intuition/bootstrap.ipynb @@ -9,7 +9,7 @@ "\n", "**Students tasks**\n", "\n", - "## Task A : Descriptive stats\n", + "## Task A : Descriptive stats (Population=500,000)\n", "\n", "- [ ] Compute and report: mean, median, SD, IQR, five-number\n", "- [ ] Plot (histogram and boxplot)\n", diff --git a/session-2-intuition/main.ipynb b/session-2-intuition/main.ipynb new file mode 100644 index 0000000..cb978e9 --- /dev/null +++ b/session-2-intuition/main.ipynb @@ -0,0 +1,248 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "799ee5dc", + "metadata": {}, + "source": [ + "# Week 2 lab: Descriptive stattistics + Bootstrap(sampling variability)\n", + "\n", + "Keypoints:\n", + "- Descriptive statistics (center, spread, shape)\n", + "- population vs sample\n", + "- sampling variability\n", + "\n", + "Lab focus: \n", + "- Bootstrap intuition (mean vs median)\n", + "- Why estimates \"move\"\n", + "- effect of sample size" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "3cd42504", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0.00123015, 0.29874554, -0.27413786])" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# cell 1: Imports + settings\n", + "import numpy as np \n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "\n", + "rng = np.random.default_rng(7) \n", + "rng.normal(size=3)" + ] + }, + { + "cell_type": "markdown", + "id": "7dee5625", + "metadata": {}, + "source": [ + "# Example 1: Descriptive statistic\n", + "We will compute the mean, median, SD, variance, IQR, five-number summary, and flag out the outliers using the $1.5\\times IQR$ rule." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "7e5313f5", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(0 2.720289\n", + " 1 3.251926\n", + " 2 2.306007\n", + " 3 1.593041\n", + " 4 2.069273\n", + " Name: x, dtype: float64,\n", + " count 60.000000\n", + " mean 3.025476\n", + " std 2.512343\n", + " min 0.600462\n", + " 25% 1.674114\n", + " 50% 2.626124\n", + " 75% 3.280787\n", + " max 18.123107\n", + " Name: x, dtype: float64)" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Cell 2: Build a dataset with skew + outlier\n", + "rng = np.random.default_rng(7) \n", + "n = 60 \n", + "x = rng.lognormal(mean=1.0, sigma=0.6, size=n) #right-skew\n", + "x[-1] *= 10\n", + "x = pd.Series(x, name=\"x\")\n", + "x.head(), x.describe()" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "34a14a88", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "({'min': np.float64(0.6004620561861885),\n", + " 'Q1': np.float64(1.6741136607158649),\n", + " 'median': np.float64(2.6261243464658732),\n", + " 'Q3': np.float64(3.2807868747809836),\n", + " 'max': np.float64(18.12310679128717)},\n", + " np.float64(1.6066732140651188),\n", + " np.float64(-0.7358961603818135),\n", + " np.float64(5.690796695878662))" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Cell 3: Five-numbers + IQR + outliers fence:\n", + "q1 = x.quantile(0.25)\n", + "q2 = x.quantile(0.50)\n", + "q3 = x.quantile(0.75)\n", + "iqr = q3-q1\n", + "\n", + "lowerFence = q1-1.5*iqr\n", + "upperFence = q3+1.5*iqr\n", + "\n", + "fiveN = {\n", + " \"min\": x.min(),\n", + " \"Q1\": q1,\n", + " \"median\": x.median(),\n", + " \"Q3\": q3,\n", + " \"max\": x.max()\n", + "}\n", + "fiveN, iqr, lowerFence, upperFence" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "id": "1886a0a7", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "({'mean': np.float64(3.0254764501575737),\n", + " 'median': np.float64(2.6261243464658732),\n", + " 'std': np.float64(2.512343173247777),\n", + " 'var': np.float64(6.31186822016471),\n", + " 'IQR': np.float64(1.6066732140651188),\n", + " 'nOutliers': 5},\n", + " 7 6.074679\n", + " 44 6.142882\n", + " 49 9.027269\n", + " 58 6.443769\n", + " 59 18.123107\n", + " Name: x, dtype: float64)" + ] + }, + "execution_count": 31, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "outliers = x[(xupperFence)]\n", + "summary = {\n", + " \"mean\": x.mean(),\n", + " \"median\": x.median(),\n", + " \"std\": x.std(ddof=1),\n", + " \"var\": x.var(ddof=1),\n", + " \"IQR\": iqr,\n", + " \"nOutliers\": len(outliers)\n", + "}\n", + "summary, outliers" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "id": "af45b496", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAHHCAYAAACle7JuAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJT5JREFUeJzt3QmUFdWBBuDLDhpoFFRAZXEBFRBU1IlmjAsBlVEYExUHDWpcxuDCmBhhDCrREVyCnCiDy4yiRlBjRDxicATBZcSI4IImoigSEhdcaRFpHHhzbp3TfWjoZuk0/d59/X3nFN1Vr6r61qu3/NylqkEul8sFAIAENcx3AQAAakqQAQCSJcgAAMkSZACAZAkyAECyBBkAIFmCDACQLEEGAEiWIAMAJEuQAWpF586dw5lnnunZBOqUIANUadKkSaFBgwbh5ZdfrvLxI488MvTo0ePvevaeeOKJcPXVVzsDQI0JMkCtWLRoUbjzzju3OsiMHj3aGQBqTJABakWzZs1CkyZNkno2v/7663wXAfg7CTLANukj8+2332a1LXvvvXdo3rx5aNOmTfje974XnnrqqezxuO6ECROy32MTVvm0fsj42c9+FnbfffcsJHXr1i3cdNNNIZfLVfq733zzTbj44otD27ZtQ8uWLcOJJ54Y/va3v2X7Wr/ZKv4el/3pT38K//Iv/xJ22GGHrDzR66+/npVnjz32yMrarl27cPbZZ4fPPvus0t8q38fbb78dTj/99FBSUhJ22mmnMGrUqKxcy5YtCwMHDgytWrXK9vHrX//aqwu2scbb+g8AaVuxYkX49NNPN1oeg8qmxC/9MWPGhHPOOScccsghobS0NOtvs2DBgvCDH/wgnH/++eGDDz7Igs19991XadsYCmIgmT17dvjJT34SevfuHZ588slw2WWXZSHl5ptvrlg3BpCHHnoonHHGGeEf/uEfwjPPPBMGDBhQbblOPvnkLFxdd911FaEoluG9994LZ511VhZA3nzzzXDHHXdkP1988cVKASs69dRTw7777hvGjh0bpk+fHq699tqw4447httvvz0cffTR4frrrw/3339/+PnPfx4OPvjgcMQRR2zx8w1spRxAFe6+++74Lb/JqXv37hXrd+rUKTd06NCK+V69euUGDBiwyed22LBh2X429Oijj2bLr7322krLf/SjH+UaNGiQW7x4cTY/f/78bL3hw4dXWu/MM8/Mll911VUVy+Lvcdlpp5220d9btWrVRsumTJmSrf/ss89utI/zzjuvYtn//d//5XbbbbesXGPHjq1Y/sUXX+RatGhR6TkBap+mJWCTYvNPrLHYcNp///03uV3r1q2zGo133nlnq5/h2Am4UaNGWZPR+mJTU6xF+cMf/pDNz5gxI/v505/+tNJ6F110UbX7/td//deNlrVo0aLi99WrV2c1ULF2J4o1SBuKtUzlYjn79OmTlSvWHq1//LE5LNb0ANuOpiVgk2KzUPyi3lDsY1JVk1O5X/3qV1l/ka5du2bDtI899tis+WdzAShaunRp6NChQ9bnZX2xOaf88fKfDRs2DF26dKm03l577VXtvjdcN/r888+z/jwPPPBAWL58+UZNaxvq2LFjpfnYVyb2rYn9dDZcvmE/G6B2qZEBtonYL+Tdd98Nd911VxZk/uu//isceOCB2c98Wr/2pdwpp5ySDR2PtTWPPPJI+J//+Z+K2p5169ZttH6shdmSZdGGnZOB2iXIANtM7AAbO9BOmTIlG9ETa2PWH0m0YSfacp06dco6An/11VeVlr/11lsVj5f/jEFjyZIlldZbvHjxFpfxiy++CLNmzQojRozIamX++Z//OeuMHEcwAYVPkAG2iQ2bVL7zne9kTT5lZWUVy7bffvvs55dffllp3eOPPz6sXbs23HrrrZWWx9FKMfwcd9xx2Xz//v2zn//5n/9Zab1bbrlli8tZXpOyYc3J+PHjt3gfQP7oIwNsE/vtt192G4ODDjooq5mJQ68ffvjhcOGFF1asEx+LYqfeGEpiqBg8eHA44YQTwlFHHRWuuOKK8P7774devXplzT3Tpk0Lw4cPD3vuuWfF9j/84Q+z0BGDU/nw63idl03V+KwvXvMlNoPdcMMN2ZDyXXfdNftbG9byAIVJkAG2iRhOHnvssSwUxFqY2AwUr7cSrwVT7qSTTspGGMVOtr/97W+zWpEYZGIH3rjtlVdeGR588MFw9913Zxfcu/HGG7ORS+u79957s2u/xOarqVOnhr59+2bbxBFDsQPulpg8eXJWjjhCK5ahX79+2cio2OEYKGwN4hjsfBcCoDa9+uqr4YADDsjC0ZAhQzy5UMT0kQGSFm9RsKHY1BRrdVxRF4qfpiUgabFvy/z587M+NY0bN86ahOJ03nnnZfdpAoqbpiUgafEqw3HYdLwZ5MqVK7OL1cUL78WOwjHYAMVNkAEAkqWPDACQLEEGAEhW0Tcgx8uXx0udx5vPbcnFsQCA/ItXh4m3KYnXc4qjEOttkIkhxsgFAEhTvE/bbrvtVn+DTKyJKX8i4qXIAYDCV1pamlVElH+P19sgU96cFEOMIAMAadlctxCdfQGAZAkyAECyBBkAIFmCDACQLEEGAEiWIAMAJEuQAQCSJcgAAMkSZACAZAkyAECyBBkAIFmCDACQLEEGAEiWIAMAJEuQAQCS1TjfBaB2dR4xvaCe0vfHDsh3EQAoYmpkAIBkCTIAQLIEGQAgWYIMAJAsQQYASJYgAwAkS5ABAJIlyAAAyRJkAIBkCTIAQLIEGQAgWYIMAJAsQQYASJYgAwAkS5ABAJIlyAAAyRJkAIBkCTIAQLIEGQAgWYIMAJAsQQYASJYgAwAkK69B5tlnnw0nnHBC6NChQ2jQoEF49NFHKx779ttvw+WXXx569uwZtt9++2ydH//4x+GDDz7IZ5EBgAKS1yDz9ddfh169eoUJEyZs9NiqVavCggULwqhRo7KfjzzySFi0aFE48cQT81JWAKDwNM7nHz/uuOOyqSolJSXhqaeeqrTs1ltvDYccckj4y1/+Ejp27FhHpQQAClVeg8zWWrFiRdYE1bp162rXKSsry6ZypaWldVQ6AKCuJRNkVq9enfWZOe2000KrVq2qXW/MmDFh9OjRISWdR0zPdxEAIElJjFqKHX9POeWUkMvlwsSJEze57siRI7Oam/Jp2bJldVZOAKBuNU4lxCxdujQ8/fTTm6yNiZo1a5ZNAEDxa5xCiHnnnXfC7NmzQ5s2bfJdJACggOQ1yKxcuTIsXry4Yn7JkiXh1VdfDTvuuGNo3759+NGPfpQNvX788cfD2rVrw0cffZStFx9v2rRpHksOAIT6HmRefvnlcNRRR1XMX3rppdnPoUOHhquvvjo89thj2Xzv3r0rbRdrZ4488sg6Li0AUGjyGmRiGIkdeKuzqccAAJIYtQQAUBVBBgBIliADACRLkAEAkiXIAADJEmQAgGQJMgBAsgQZACBZggwAkCxBBgBIliADACRLkAEAkiXIAADJEmQAgGQJMgBAsgQZACBZggwAkCxBBgBIliADACRLkAEAkiXIAADJEmQAgGQJMgBAsgQZACBZggwAkCxBBgBIliADACRLkAEAkiXIAADJEmQAgGQJMgBAsgQZACBZggwAkCxBBgBIliADACRLkAEAkiXIAADJEmQAgGQJMgBAsgQZACBZggwAkCxBBgBIliADACQrr0Hm2WefDSeccELo0KFDaNCgQXj00UcrPZ7L5cKVV14Z2rdvH1q0aBH69u0b3nnnnbyVFwAoLHkNMl9//XXo1atXmDBhQpWP33DDDeE3v/lNuO2228If//jHsP3224f+/fuH1atX13lZAYDC0ziff/y4447LpqrE2pjx48eHX/7yl2HgwIHZsnvvvTfssssuWc3N4MGD67i0AEChKdg+MkuWLAkfffRR1pxUrqSkJBx66KFh7ty51W5XVlYWSktLK00AQHEq2CATQ0wUa2DWF+fLH6vKmDFjssBTPu2+++7bvKwAQH4UbJCpqZEjR4YVK1ZUTMuWLct3kQCA+hZk2rVrl/38+OOPKy2P8+WPVaVZs2ahVatWlSYAoDgVbJDp0qVLFlhmzZpVsSz2d4mjl7773e/mtWwAQGHI66illStXhsWLF1fq4Pvqq6+GHXfcMXTs2DEMHz48XHvttWHvvffOgs2oUaOya84MGjQon8UGAApEXoPMyy+/HI466qiK+UsvvTT7OXTo0DBp0qTwi1/8IrvWzHnnnRe+/PLL8L3vfS/MmDEjNG/ePI+lBgAKRYNcvGBLEYvNUXH0Uuz4W6j9ZTqPmB6K1ftjB+S7CAAU8fd3wfaRAQDYHEEGAEiWIAMAJEuQAQCSJcgAAMkSZACAZAkyAECyBBkAIFmCDACQLEEGAEiWIAMAJEuQAQCSJcgAAMkSZACAZAkyAECyBBkAIFmCDACQLEEGAEiWIAMAJEuQAQCSJcgAAMkSZACAZAkyAECyBBkAIFmCDACQLEEGAEiWIAMAJEuQAQCSJcgAAMkSZACAZAkyAECyBBkAIFmCDACQLEEGAEiWIAMAJEuQAQCSJcgAAMkSZACAZAkyAECyBBkAIFmCDACQLEEGAEiWIAMAJKugg8zatWvDqFGjQpcuXUKLFi3CnnvuGa655pqQy+XyXTQAoAA0DgXs+uuvDxMnTgz33HNP6N69e3j55ZfDWWedFUpKSsLFF1+c7+IBAHlW0EHmhRdeCAMHDgwDBgzI5jt37hymTJkSXnrppXwXDQAoAAXdtHTYYYeFWbNmhbfffjubf+2118Lzzz8fjjvuuGq3KSsrC6WlpZUmAKA4FXSNzIgRI7Igss8++4RGjRplfWb+4z/+IwwZMqTabcaMGRNGjx5dJ+XrPGJ6nfwdACDBGpmHHnoo3H///WHy5MlhwYIFWV+Zm266KftZnZEjR4YVK1ZUTMuWLavTMgMAdaega2Quu+yyrFZm8ODB2XzPnj3D0qVLs1qXoUOHVrlNs2bNsgkAKH4FXSOzatWq0LBh5SLGJqZ169blrUwAQOEo6BqZE044IesT07Fjx2z49SuvvBLGjRsXzj777HwXDQAoAAUdZG655Zbsgng//elPw/Lly0OHDh3C+eefH6688sp8Fw0AKAAFHWRatmwZxo8fn00AAEn1kQEA2BRBBgBIliADANSvIBMvTrdw4cKK+WnTpoVBgwaFf//3fw9r1qypzfIBANRukIkjh8rvf/Tee+9lF6zbbrvtwu9+97vwi1/8oia7BAComyATQ0zv3r2z32N4OeKII7LbCEyaNCn8/ve/r8kuAQDqJsjkcrmKq+vOnDkzHH/88dnvu+++e/j0009rsksAgLoJMn369AnXXnttuO+++8IzzzwTBgwYkC1fsmRJ2GWXXWqySwCAugkyN998c9bh98ILLwxXXHFF2GuvvbLlDz/8cDjssMNqsksAgLq5sm+vXr0qjVoqd+ONN4bGjQv6YsEAQH2vkdljjz3CZ599ttHy1atXh65du9ZGuQAAtk2Qef/998PatWs3Wl5WVhb++te/1mSXAABbbavagR577LGK35988slQUlJSMR+DzaxZs0KXLl2cBgCg8IJMvHpv1KBBgzB06NBKjzVp0iR07tw5/PrXv67dEgIA1EaQKb92TKx1mTdvXmjbtu3WbA4AUKtqNMQoXi8GACDfajxWOvaHidPy5csramrK3XXXXbVRNgCA2g8yo0ePDr/61a+yK/y2b98+6zMDAJBEkLntttuyG0SeccYZtV8iAIBteR2ZNWvWuBUBAJBmkDnnnHPC5MmTa780AADbumkp3orgjjvuCDNnzgz7779/dg2Z9Y0bN64muwUA2PZB5vXXXw+9e/fOfn/jjTcqPabjLwBQ0EFm9uzZtV8SAIC66CMDAJBsjcxRRx21ySakp59++u8pEwDAtgsy5f1jyn377bfh1VdfzfrLbHgzSQCAggoyN998c5XLr7766rBy5cq/t0wAAHXfR+b00093nyUAIM0gM3fu3NC8efPa3CUAQO02LZ100kmV5nO5XPjwww/Dyy+/HEaNGlWTXQIA1E2QKSkpqTTfsGHD0K1bt+yO2P369avJLgEA6ibI3H333TXZDAAg/0Gm3Pz588Of//zn7Pfu3buHAw44oLbKBQCwbYLM8uXLw+DBg8OcOXNC69ats2VffvlldqG8Bx54IOy000412S0AwLYftXTRRReFr776Krz55pvh888/z6Z4MbzS0tJw8cUX12SXAAB1UyMzY8aMMHPmzLDvvvtWLNtvv/3ChAkTdPYFAAq7RmbdunWhSZMmGy2Py+JjAAAFG2SOPvrocMkll4QPPvigYtnf/va38G//9m/hmGOOqc3yAQDUbpC59dZbs/4wnTt3DnvuuWc2denSJVt2yy231GSXAAB100dm9913DwsWLMj6ybz11lvZsthfpm/fvjXZHQDAtq+Refrpp7NOvbHmpUGDBuEHP/hBNoIpTgcffHB2LZnnnnuuZiUBANiWQWb8+PHh3HPPDa1atarytgXnn39+GDdu3NaWAQBg2weZ1157LRx77LHVPh7vsxSv9lubYifi008/PbRp0ya0aNEi9OzZM7s5JQDAVvWR+fjjj6scdl2ucePG4ZNPPqm1Z/WLL74Ihx9+eHbF4D/84Q/ZFYPfeeedsMMOO9Ta3wAA6kmQ2XXXXbMr+O61115VPv7666+H9u3b11bZwvXXX591LF7/JpVxdBQAwFY3LR1//PFh1KhRYfXq1Rs99s0334Srrroq/NM//VOtPbOPPfZY6NOnTzj55JPDzjvvnN2U8s4779zkNmVlZVln5PUnAKA4bVWQ+eUvf5ndV6lr167hhhtuCNOmTcumWHPSrVu37LErrrii1gr33nvvhYkTJ4a99947PPnkk+GCCy7I7uV0zz33VLvNmDFjso7H5VOs0QEAilODXC6X25oNli5dmgWKGCzKN41Dsfv375/da6k2m36aNm2a1ci88MILFctikJk3b16YO3dutTUycSoXa2RimFmxYkWVo63+Hp1HTK/V/RWj98cOyHcRAEhQ/P6OFRKb+/7e6gviderUKTzxxBNZR9zFixdnYSbWmGyLDrixv028bs364oX3fv/731e7TbNmzbIJACh+NbqybxSDS7wI3rYURywtWrSo0rK33347C1MAADW611JdiTehfPHFF8N1112X1f5Mnjw53HHHHWHYsGH5LhoAUAAKOsjEGp+pU6eGKVOmhB49eoRrrrkmu7rwkCFD8l00ACDlpqW6Eodz1+aQbgCgeBR0jQwAwKYIMgBAsgQZACBZggwAkCxBBgBIliADACRLkAEAkiXIAADJEmQAgGQJMgBAsgQZACBZggwAkCxBBgBIliADACRLkAEAkiXIAADJEmQAgGQJMgBAsgQZACBZggwAkCxBBgBIliADACRLkAEAkiXIAADJEmQAgGQJMgBAsgQZACBZggwAkCxBBgBIliADACRLkAEAkiXIAADJEmQAgGQJMgBAsgQZACBZggwAkCxBBgBIliADACRLkAEAkiXIAADJEmQAgGQJMgBAsgQZACBZSQWZsWPHhgYNGoThw4fnuygAQAFIJsjMmzcv3H777WH//ffPd1EAgAKRRJBZuXJlGDJkSLjzzjvDDjvskO/iAAAFIokgM2zYsDBgwIDQt2/fza5bVlYWSktLK00AQHFqHArcAw88EBYsWJA1LW2JMWPGhNGjR2/zcrFlOo+YXlBP1ftjB+S7CADUlxqZZcuWhUsuuSTcf//9oXnz5lu0zciRI8OKFSsqprgPAKA4FXSNzPz588Py5cvDgQceWLFs7dq14dlnnw233npr1ozUqFGjSts0a9YsmwCA4lfQQeaYY44JCxcurLTsrLPOCvvss0+4/PLLNwoxAED9UtBBpmXLlqFHjx6Vlm2//fahTZs2Gy0HAOqfgu4jAwCQbI1MVebMmZPvIgAABUKNDACQLEEGAEiWIAMAJEuQAQCSJcgAAMkSZACAZAkyAECyBBkAIFmCDACQLEEGAEiWIAMAJEuQAQCSJcgAAMkSZACAZAkyAECyBBkAIFmCDACQLEEGAEiWIAMAJEuQAQCSJcgAAMkSZACAZAkyAECyBBkAIFmCDACQLEEGAEiWIAMAJEuQAQCSJcgAAMkSZACAZAkyAECyBBkAIFmCDACQLEEGAEiWIAMAJEuQAQCSJcgAAMkSZACAZAkyAECyBBkAIFmCDACQLEEGAEiWIAMAJKugg8yYMWPCwQcfHFq2bBl23nnnMGjQoLBo0aJ8FwsAKBAFHWSeeeaZMGzYsPDiiy+Gp556Knz77behX79+4euvv8530QCAAtA4FLAZM2ZUmp80aVJWMzN//vxwxBFH5K1cAEBhKOgamQ2tWLEi+7njjjvmuygAQAEo6BqZ9a1bty4MHz48HH744aFHjx7VrldWVpZN5UpLS+uohABAXUsmyMS+Mm+88UZ4/vnnN9tBePTo0XVWLuqnziOm18p+3h87oFb2A1BfJdG0dOGFF4bHH388zJ49O+y2226bXHfkyJFZE1T5tGzZsjorJwBQtwq6RiaXy4WLLrooTJ06NcyZMyd06dJls9s0a9YsmwCA4te40JuTJk+eHKZNm5ZdS+ajjz7KlpeUlIQWLVrku3gAQJ4VdNPSxIkTs+ahI488MrRv375ievDBB/NdNACgABR80xIAQJI1MgAAmyLIAADJEmQAgGQJMgBAsgQZACBZggwAkCxBBgBIliADACRLkAEAkiXIAADJEmQAgGQJMgBAsgQZACBZggwAkCxBBgBIliADACRLkAEAkiXIAADJEmQAgGQJMgBAsgQZACBZjfNdAKhLnUdML8ryvD92QChGnh/w/tocNTIAQLIEGQAgWYIMAJAsQQYASJYgAwAkS5ABAJIlyAAAyRJkAIBkCTIAQLIEGQAgWYIMAJAsQQYASJYgAwAkS5ABAJIlyAAAyRJkAIBkCTIAQLIEGQAgWYIMAJAsQQYASJYgAwAkK4kgM2HChNC5c+fQvHnzcOihh4aXXnop30UCAApAwQeZBx98MFx66aXhqquuCgsWLAi9evUK/fv3D8uXL8930QCAPCv4IDNu3Lhw7rnnhrPOOivst99+4bbbbgvbbbdduOuuu/JdNAAgzwo6yKxZsybMnz8/9O3bt2JZw4YNs/m5c+fmtWwAQP41DgXs008/DWvXrg277LJLpeVx/q233qpym7Kysmwqt2LFiuxnaWlprZdvXdmqWt8n1MS2eH0Xgtp6jxXr8wPF/P4q328ul0s3yNTEmDFjwujRozdavvvuu+elPFAXSsZ7nj0/UJyfP1999VUoKSlJM8i0bds2NGrUKHz88ceVlsf5du3aVbnNyJEjs87B5datWxc+//zz0KZNm9CgQYON0l4MOMuWLQutWrUK9U19P/6ovj8Hjr9+n//Ia6B+vwZKC/j4Y01MDDEdOnTY5HoFHWSaNm0aDjrooDBr1qwwaNCgimAS5y+88MIqt2nWrFk2ra9169ab/Dvx5BXaCaxL9f34o/r+HDj++n3+I6+B+v0aaFWgx7+pmpgkgkwUa1eGDh0a+vTpEw455JAwfvz48PXXX2ejmACA+q3gg8ypp54aPvnkk3DllVeGjz76KPTu3TvMmDFjow7AAED9U/BBJorNSNU1Jf09YhNUvNDehk1R9UV9P/6ovj8Hjr9+n//Ia6B+vwaK4fw3yG1uXBMAQIEq6AviAQBsiiADACRLkAEAkiXIAADJKvogM2HChNC5c+fQvHnzcOihh4aXXnppk+v/7ne/C/vss0+2fs+ePcMTTzwRUr1Vw8EHHxxatmwZdt555+yCgosWLdrkNpMmTcqufrz+FJ+HVF199dUbHU88t/Xh/Efxdb/h8cdp2LBhRXv+n3322XDCCSdkVwKN5X/00UcrPR7HNsRLObRv3z60aNEiuwHtO++8U+ufI4V4/N9++224/PLLs9f19ttvn63z4x//OHzwwQe1/j4q1PN/5plnbnQsxx57bNGc/y15Dqr6TIjTjTfeGFJ9DRR1kHnwwQezC+rFoWULFiwIvXr1Cv379w/Lly+vcv0XXnghnHbaaeEnP/lJeOWVV7Iv/zi98cYbITXPPPNM9oX14osvhqeeeir7EOvXr192McFNiVd2/PDDDyumpUuXhpR179690vE8//zz1a5bTOc/mjdvXqVjj6+D6OSTTy7a8x9f3/F9Hr94qnLDDTeE3/zmN+G2224Lf/zjH7Mv9PiZsHr16lr7HCnU41+1alVW/lGjRmU/H3nkkew/NyeeeGKtvo8K+fxHMbisfyxTpkzZ5D5TOv9b8hysf+xxuuuuu7Jg8sMf/jAk+xrIFbFDDjkkN2zYsIr5tWvX5jp06JAbM2ZMleufcsopuQEDBlRaduihh+bOP//8XOqWL18eh9nnnnnmmWrXufvuu3MlJSW5YnHVVVflevXqtcXrF/P5jy655JLcnnvumVu3bl29OP/x9T516tSK+Xjc7dq1y914440Vy7788stcs2bNclOmTKm1z5FCPf6qvPTSS9l6S5curbX3USEf/9ChQ3MDBw7cqv2kev639DUQn4+jjz56k+sU+mugaGtk1qxZE+bPn59VHZdr2LBhNj937twqt4nL118/ism7uvVTsmLFiuznjjvuuMn1Vq5cGTp16pTdRGzgwIHhzTffDCmLzQaxinWPPfYIQ4YMCX/5y1+qXbeYz398P/z2t78NZ5999kY3Ty3m87++JUuWZFcHX/8cx/u4xKaC6s5xTT5HUvtciK+Hzd2PbmveR4Vuzpw5WXN7t27dwgUXXBA+++yzatct9vP/8ccfh+nTp2e10JtTyK+Bog0yn376aVi7du1GtzKI8/HDrCpx+dasn4p4o83hw4eHww8/PPTo0aPa9eIbO1YzTps2LfvSi9sddthh4a9//WtIUfyCiv0+4i0tJk6cmH2R/eM//mN2N9X6dP6j2E7+5ZdfZn0E6sv531D5edyac1yTz5FUxOa02GcmNqdu6maBW/s+KmSxWenee+/Nbjx8/fXXZ03wxx13XHaO69v5j+65556sH+VJJ50UNqXQXwNJ3KKAv0/sKxP7eWyuTfO73/1uNpWLX2L77rtvuP3228M111yT3GmIH1Dl9t9//+zNGGsbHnrooS36H0gx+e///u/s+Yj/o6ov55/qxT5zp5xyStb5OX4x1Zf30eDBgyt+j52e4/HsueeeWS3NMcccE+qbu+66K6td2Vyn/kJ/DRRtjUzbtm1Do0aNsqqz9cX5du3aVblNXL4166cg3qPq8ccfD7Nnzw677bbbVm3bpEmTcMABB4TFixeHYhCrz7t27Vrt8RTj+Y9ih92ZM2eGc845p16f//LzuDXnuCafI6mEmPi6iB3AN1UbU5P3UUpiM0k8x9UdSzGe/3LPPfdc1tl7az8XCvE1ULRBpmnTpuGggw7KqhDLxaryOL/+/zrXF5evv34U3+jVrV/I4v+0YoiZOnVqePrpp0OXLl22eh+xSnXhwoXZUNViEPt/vPvuu9UeTzGd//XdfffdWZ+AAQMG1OvzH98D8ctn/XNcWlqajV6q7hzX5HMkhRAT+zvEcNumTZtafx+lJDabxj4y1R1LsZ3/DWtp47HFEU7JvwZyReyBBx7IRiRMmjQp96c//Sl33nnn5Vq3bp376KOPssfPOOOM3IgRIyrW/9///d9c48aNczfddFPuz3/+c9ZTu0mTJrmFCxfmUnPBBRdkI1DmzJmT+/DDDyumVatWVayz4fGPHj069+STT+befffd3Pz583ODBw/ONW/ePPfmm2/mUvSzn/0sO/4lS5Zk57Zv3765tm3bZiO4iv38rz/ComPHjrnLL798o8eK8fx/9dVXuVdeeSWb4sfbuHHjst/LR+WMHTs2+wyYNm1a7vXXX89GbHTp0iX3zTffVOwjjuC45ZZbtvhzJJXjX7NmTe7EE0/M7bbbbrlXX3210udCWVlZtce/ufdRKscfH/v5z3+emzt3bnYsM2fOzB144IG5vffeO7d69eqiOP9b8h6IVqxYkdtuu+1yEydOzFUltddAUQeZKJ6M+EHetGnTbBjdiy++WPHY97///Ww43voeeuihXNeuXbP1u3fvnps+fXouRfEFXNUUh9hWd/zDhw+veK522WWX3PHHH59bsGBBLlWnnnpqrn379tnx7Lrrrtn84sWL68X5LxeDSTzvixYt2uixYjz/s2fPrvJ1X36ccQj2qFGjsuOLX07HHHPMRs9Np06dshC7pZ8jqRx//BKq7nMhblfd8W/ufZTK8cf/xPXr1y+30047Zf9Bicd57rnnbhRIUj7/W/IeiG6//fZcixYtsssPVCW110CD+E++a4UAAGqiaPvIAADFT5ABAJIlyAAAyRJkAIBkCTIAQLIEGQAgWYIMAJAsQQYASJYgAwAkS5ABAJIlyABJ+eSTT7K7WF933XUVy1544YXsTsUb3r0cKH7utQQk54knngiDBg3KAky3bt1C7969w8CBA8O4cePyXTSgjgkyQJKGDRsWZs6cGfr06RMWLlwY5s2bF5o1a5bvYgF1TJABkvTNN9+EHj16hGXLloX58+eHnj175rtIQB7oIwMk6d133w0ffPBBWLduXXj//ffzXRwgT9TIAMlZs2ZNOOSQQ7K+MbGPzPjx47PmpZ133jnfRQPqmCADJOeyyy4LDz/8cHjttdfCd77znfD9738/lJSUhMcffzzfRQPqmKYlIClz5szJamDuu+++0KpVq9CwYcPs9+eeey5MnDgx38UD6pgaGQAgWWpkAIBkCTIAQLIEGQAgWYIMAJAsQQYASJYgAwAkS5ABAJIlyAAAyRJkAIBkCTIAQLIEGQAgWYIMABBS9f8wCSwmeyOxIAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAHHCAYAAADjzRHEAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAHP5JREFUeJzt3QnwVlX9P/DDIpsCjrmBG2igjAvulqbmMrggy5iamoWlZoaVS26VouVoLpktjpmTYDJZ6riMipqYiOYauIsLhGaBipaAAiJw/3Pub57v/8su8OHLs7xeMw+X5977LPecu7y/5557n1ZFURQJACBA64g3AQAQLACAUFosAIAwggUAEEawAADCCBYAQBjBAgAII1gAAGEECwAgjGABrDGtWrVKF154oRqAOiJYQJ0aMWJEeeBu/thwww3Tfvvtl+67775Uy1555ZUykLz55ptr+qsAi2i76Aigvvz0pz9NPXv2TPlngd59990ycBx66KHp7rvvTocddliq1WBx0UUXpS9/+cupR48ea/rrAM0IFlDnDjnkkLTrrrs2PT/hhBPSRhttlG6++eaaDRZA9XIqBBrMuuuumzp27Jjatv3/f1d8/PHH6cwzz0ybbbZZat++fdp6663TlVdeWbZyZLNnz07bbLNN+cj/r/jvf/+bunXrlvbcc880f/78ctzxxx+f1llnnfTPf/4zHXTQQWnttddO3bt3L1tOPsuPKT/77LNlGOrSpUv5PgcccEB68sknm6bnFpcjjzyy/H8+rVM5zTNmzJjQcgJWjmABdW769Onp/fffT9OmTUsvv/xyOuWUU9JHH32UjjvuuHJ6PtgPHDgw/fKXv0wHH3xwuuqqq8pgcdZZZ6UzzjijnCcHkRtvvDFNnDgx/fjHP25676FDh5bvnw/2bdq0aRqfQ0Z+r9wycvnll6dddtklDRs2rHwsS/5+e++9d3r++efT2Wefnc4///w0efLk8pTHU089Vc6zzz77pO9///vl/3/0ox+lm266qXz06dNntZQfsIIKoC4NHz48Nw8s9mjfvn0xYsSIpvnuvPPOcvzFF1+80OuPOOKIolWrVsXEiRObxp133nlF69ati7Fjxxa33npr+bqrr756odcNGTKkHP+9732vadyCBQuK/v37F+3atSumTZvWND7PN2zYsKbngwcPLueZNGlS07gpU6YUnTt3LvbZZ5+mcZXPfvjhh0PKCoijxQLq3DXXXJMefPDB8jFy5Mjy9MGJJ56Ybr/99nL6qFGjytaGSitART41ko/9za8gyVdibLvttmnIkCHpu9/9btp3330Xe13Fqaee2vT/fKoiP587d24aPXr0EufPrRx//etf0+DBg9OWW27ZND6fajn22GPTY489lmbMmLHK5QGsXjpvQp3bfffdF+q8ecwxx6SddtqpPNDnzptvvfVW2Qeic+fOC72ucmohT69o165duuGGG9Juu+2WOnTokIYPH16GhkW1bt16oXCQ9e7duxwu7RLRfKpm1qxZ5WmYReXvsmDBgvT222+XwQaoXlosoMHkg35utZg6dWp64403Vvj1DzzwQDmcM2fOSr0eqG+CBTSgefPmlcPciXOLLbZIU6ZMSTNnzlxonldffbUc5ukVL7zwQnl1xze/+c2y1SOfUsmdNxeVWxfyVSHNvf766+Vwafed2GCDDVKnTp3Sa6+9tti0/F1yIMpXrWRLaiUBqoNgAQ3m008/Lfsy5NMa+RRDvllW7t/w29/+dqH58lUi+QCeL/2svC5fSppPm/zqV78qrwTJN9w6/fTTl/g5zd8v99XIz9daa63y8tElyf08+vXrl+66666FTpfkz/jTn/6UvvSlL5WXoGb5Etbsww8/DCgRIJI+FlDncufLSuvDe++9Vx6k8ymMc889tzxQDxgwoDw1ki8jzQf0vn37lsEjH+BPO+20tNVWW5Wvvfjii9Nzzz2XHnroobI/xg477JAuuOCC9JOf/CQdccQRZUCpyP0v7r///rKT5x577FF+h3vvvbe8PDS3TCxN/ozcyTSHiNw5NN9r47rrrkuffPJJedlqxY477lgGkcsuu6xsMcn33th///3LW5YDa1jgFSZAlV9u2qFDh2LHHXcsrr322vIS0IqZM2cWp59+etG9e/dirbXWKnr16lVcccUVTfOMGzeuaNu27UKXkGbz5s0rdtttt/J1//vf/5ouN1177bXLS0b79etXdOrUqdhoo43Ky0rnz5+/0OsXvdw0Gz9+fHHQQQcV66yzTvna/fbbr3j88ccXW77rr7++2HLLLYs2bdq49BSqSKv8z5oON0D9yKdLbrvttrL/BtB49LEAAMIIFgBAGMECAAijjwUAEEaLBQAQRrAAAGr3Bln5Vr/59sH5BjtuywsAtSHfnSLf+j/ffTffYr9qgkUOFZX7/QMAtSX/yvCmm25aPcGi8tPM+YtV7vsPAFS3GTNmlA0DleN41QSLyumPHCoECwCoLcvrxqDzJgAQRrAAAMIIFgBAGMECAAgjWAAAYQQLACCMYAEAhBEsAIAwggUAEEawAADCCBYAgGABAFQfLRYAQBjBAgAII1gAAGEECwAgjGABAIQRLACAMIIFABBGsAAAwggWAEAYwQIACCNYAABhBAsAIIxgAQCEESwAgDCCBQAQRrAAAMIIFgBAGMECAAgjWAAAYQQLACCMYAEAhBEsAIAwggUAEEawAADCCBYAQBjBAgAII1gAAGEECwAgjGABAIQRLACAMIIFABBGsAAAwggWAEAYwQIACCNYAABhBAsAIIxgAQCEESwAgDCCBQAQRrAAAMIIFgBAGMECAAgjWAAAYQQLACCMYAEACBYAQPXRYgEAhBEsAIAwggUAEEawAADCCBYAQBjBAgAII1gAAGEECwAgjGABAIQRLACAMIIFABBGsAAAwggWAEAYwQIACCNYAABhBAsAIIxgAQCEESwAgDCCBQAQRrAAAMIIFgBAGMECAAgjWAAAYQQLACCMYAEAhBEsAIAwggUAEEawAADCCBYAQBjBAgAII1gAAGEECwAgjGABAIQRLACAMIIFABBGsAAAwggWAEAYwQIACCNYAABhBAsAIIxgAQCEESwAgDCCBQAQRrAAAMIIFgCAYAEAVB8tFgBAGMECAAgjWAAAYQQLACCMYAEAhBEsAIAwggUAEEawAADCCBYAQBjBAgAII1gAAGEECwAgjGABAIQRLACAMIIFABBGsAAAwggWAEAYwQIACCNYAABhBAsAIIxgAQCEESwAgDCCBQAQRrAAAMIIFgBAGMECAAgjWAAAYQQLACCMYAEAhBEsAIAwggUAEEawAADCCBYAQBjBAgAII1gAAGEECwAgjGABAIQRLACAMIIFABBGsAAAwggWAEAYwQIACCNYAABhBAsAIIxgAQAIFgBA9dFiAQCEESwAgDCCBQAQpm3cWzWeN954I82cOTNVs86dO6devXqt6a8BQIMQLFYhVPTu3TutLhuv0yqdvEu7dN24uemdj4pVeq/XX39duACgRQgWK6nSUjFy5MjUp0+fFK3jh6+nPmNPTl+9YESave7KBZgJEyak4447rupbVQCoH4LFKsqhYuedd07hprROaWxKfbbZJqXuO8a/PwCsBjpvAgBhBAsAIIxgAQCEESwAgDCCBQAQRrAAAMIIFgBAGMECAAgjWAAAYQQLACBM3QSLWbNmpfHjx5dD6oM6Bag9dRMsXn311bTLLruUQ+qDOgWoPXUTLACANU+wAADCCBYAQBjBAgAII1gAAGEECwAgjGABAIQRLACAMG3j3gpqx/z589Ojjz6a/vOf/6Rp06alDTbYIG2yySZp7733Tm3atFlsvqlTp6Zu3botNr1RKAeofvOrZH+1wi0WY8eOTQMGDEjdu3dPrVq1Snfeeefq+WZ8Zk9MeSINunNQOWT5br/99vT5z38+7bfffum4445Lp59+ejnMz/P4PH3R+Y499tjFpjcK5QDV7/Yq2l+tcLD4+OOPU9++fdM111yzer4RK6QoivSr8b9K/5z+z3KYn7N0eSM74ogj0vrrr18G40MOOSRdf/315TDL4/P0s88+uxxuv/326YknnkgzZ84sh/l5Ht8o4aJSXo1eDlDNbq+27bRYBfnld9xxxwq9Zvr06eXr8jDSuHHjyvfNw5aw2j/vP88WxbAu/zdchsf+/Vix3Yjtmh75eYt9x9Us+vvPmzev6NGjR3HYYYeVwwEDBhTz588vp+Vhfp7H9+/fv2jbtm05X2V6RWW+nj17lu9Xzyrl1bycGrEcoJrNa8Ht9LMev1d7H4tPPvmkfFTMmDFjtXzO7Nmzy+GECRNSS6h8TuVz14Sc7X7z7G9S61at04JiQTnMz/fsvmf513hLl0m1l3E+9/jmm2+ms846K91zzz3p5ptvTq1b/1+jXR6ed955ac8990wDBw5M8+bNK1sxKtMrms+X3+/LX/5yqleV8mpeTo1YDlDNHq3C7XS1B4tLL700XXTRRav7Y8qCzfK58paUP3evvfZKa8LjUx5PL3/wctPzHC7y8zx+r032WmNlUq1lnDs0ZR07diyH22233ULTK8/nzJmz0HyLqsxXeb96VVm+Rcup0coBqtnUKtxOV3uwyGnpjDPOWKjFYrPNNgv/nB49epTDkSNHpj59+qSW+Gs6H7Arn7umWysqmrdatHSZVHsZ517SzVtAXnrppfSFL3yhaXp+nnXo0GGh+RZVma/yfvWqsnyLllOjlQNUs25VuJ2u9mDRvn378rG6Vf66zAfQnXfeebV/3qKfu6ZbK5bUarGmyqRayzhfepVDyn333VcOL7nkkvKqptxcuGDBgrJ1LY+fNGlSatu2bTnfd77znYWaFyvz9ezZs3y/elYpr+bl1IjlANVs7yrcTt0gqwZVWitapVZLnJ7H5+muEFlYvp77F7/4Rbr33nvLqz9yP4vDDjss/f73vy+Hd999dzl+1KhR5SWoeb7Bgwcv1Ms6P8+vu/LKK+v+fhaV8srL28jlANWsTRVupyvcYvHRRx+liRMnNj2fPHlyeu6559J6662XNt988+jvxxJ8uuDT9M7H76QiLfnS0jw+T59XzFN+izj88MPTbbfdls4888wyeOVWifyo+OCDD8rpeb7crJjnyx2fKnLyr0xvtPJq5HKAanZ4lW2nKxws/vGPf5Q33qio9J8YMmRIGjFiROy3Y4natWmX/nzYn9N/5/x3qSW0Xof10pTXpijBJcgb2aBBg5Z7583m863pO9mtScoBqt/hVbS/WuFgkS9X0cS+5m289sblY1mmJMFiafLG9lkuvfqs89U75QDVr02V7K/0sQAAwggWAEAYwQIACCNYAABhBAsAIIxgAQAIFgBA9dFiAQCEqZtgsc0226Rx48aVQ+qDOgWoPav9101bSqdOnWr6FzxZnDoFqD1102IBAKx5ggUAEEawAADCCBYAQBjBAgAII1gAAGEECwAgjGABAIQRLACAMIIFABCmbm7p3dJmzZpVDsePH79a3r/jh6+nPimlCa++mma/s2Cl3mPChAnh3wsAlkWwWEmvvvpqOTzppJPS6rDxOq3Sybu0S9f94tj0zkfFKr1X586dw74XACyLYLGSBg8e3PQLnPnHslaXgav4+hwqevXqFfRtAGDZWhVFsWp/Dq+gGTNmpK5du6bp06enLl26tORHAwCr+fit8yYAEEawAADCCBYAQBjBAgAII1gAAGEECwAgjGABAIQRLACAMIIFABBGsAAAwggWAEAYwQIACCNYAABhBAsAIIxgAQCEESwAgDCCBQAQRrAAAMIIFgBAGMECAAgjWAAAYQQLACCMYAEAhBEsAIAwggUAEEawAADCCBYAQBjBAgAQLACA6qPFAgAII1gAAGEECwAgjGABAIQRLACAMIIFABBGsAAAwggWAEAYwQIACCNYAABhBAsAIIxgAQCEESwAgDCCBQAQRrAAAMIIFgBAGMECAAgjWAAAYQQLACCMYAEAhBEsAIAwggUAEEawAADCCBYAQBjBAgAII1gAAGEECwAgjGABAIQRLACAMIIFABBGsAAAwggWAEAYwQIACCNYAABhBAsAIIxgAQCEESwAgDCCBQAQRrAAAMIIFgBAGMECAAgjWAAAYQQLACCMYAEAhBEsAADBAgCoPlosAIAwggUAEEawAADCCBYAQBjBAgAII1gAAGEECwAgjGABAIQRLACAMIIFABBGsAAAwggWAEAYwQIACCNYAABhBAsAIIxgAQCEESwAgDCCBQAQRrAAAMIIFgBAGMECAAgjWAAAYQQLACCMYAEAhBEsAIAwggUAEEawAADCCBYAQBjBAgAII1gAAGEECwAgjGABAIQRLACAMIIFABBGsAAAwggWAEAYwQIACCNYAABhBAsAIIxgAQCEESwAgDCCBQAQRrAAAMIIFgBAGMECABAsAIDqo8UCAAgjWAAAYQQLACCMYAEAhBEsAIAwggUAEEawAADCCBYAQBjBAgAII1gAAGEECwAgjGABAIQRLACAMIIFABBGsAAAwggWAEAYwQIACCNYAABhBAsAIIxgAQCEESwAgDCCBQAQRrAAAMIIFgBAGMECAAgjWAAAYQQLACCMYAEAhBEsAIAwggUAEEawAADCCBYAQJi2qYUVRVEOZ8yY0dIfDQCspMpxu3Icr5pgMXPmzHK42WabtfRHAwABx/GuXbsudXqrYnnRI9iCBQvSlClTUufOnVOrVq0WS0M5cLz99tupS5cuqdFY/sau/8w60NjrQKPXf9boZTCjipc/x4UcKrp3755at25dPS0W+ctsuummy5wnF2a1FWhLsvyNXf+ZdaCx14FGr/+s0cugS5Uu/7JaKip03gQAwggWAEB9Bov27dunYcOGlcNGZPkbu/4z60BjrwONXv9Zo5dB+zpY/hbvvAkA1K+qarEAAGqbYAEAhBEsAIAwggUAULvB4pprrkk9evRIHTp0SHvssUd6+umnlzn/rbfemrbZZpty/u233z6NGjUq1aJLL7007bbbbuUdRzfccMM0ePDg9Nprry3zNSNGjCjvTtr8kcuhFl144YWLLUuu10ao+4q83i9aBvkxdOjQuqz/sWPHpgEDBpR36cvf/c4771xoeu43fsEFF6Ru3bqljh07pgMPPDC98cYb4fuQai2DTz/9NJ1zzjnlur322muX83zjG98o70wcvS1V6zpw/PHHL7YsBx98cN2sA2OXs/xL2h/kxxVXXFHT9d+iweIvf/lLOuOMM8pLacaPH5/69u2bDjrooPTee+8tcf7HH388HXPMMemEE05Izz77bHkwzo+XXnop1ZpHHnmkPIA8+eST6cEHHyx3Kv369Usff/zxMl+X77w2derUpsdbb72VatW222670LI89thjS523nuq+4plnnllo+fN6kB155JF1Wf953c7beD4ILMnll1+efv3rX6ff/e536amnnioPrnl/MGfOnLB9SDWXwaxZs8plOP/888vh7bffXv6xMXDgwNBtqZrXgSwHiebLcvPNNy/zPWtpHfh4OcvffLnz44YbbiiDwle+8pXarv+iBe2+++7F0KFDm57Pnz+/6N69e3HppZcucf6jjjqq6N+//0Lj9thjj+Lkk08uat17772XL/MtHnnkkaXOM3z48KJr165FPRg2bFjRt2/fzzx/Pdd9xQ9+8INiq622KhYsWFD39Z/X9TvuuKPpeV7mjTfeuLjiiiuaxn344YdF+/bti5tvvjlsH1LNZbAkTz/9dDnfW2+9FbYtVfPyDxkypBg0aNAKvU+trgPpM9R/Lov9999/mfPUQv23WIvF3Llz07hx48rmzua/G5KfP/HEE0t8TR7ffP4sJ9OlzV9Lpk+fXg7XW2+9Zc730UcfpS222KL8UZpBgwall19+OdWq3MydmwS33HLL9LWvfS3961//Wuq89Vz3le1h5MiR6Vvf+tZiP8ZXr/Xf3OTJk9M777yzUB3n3yDIzdpLq+OV2YfU4n4hrw/rrrtu2LZU7caMGVOeHt56663TKaeckj744IOlzlvP68C7776b7r333rKVdnmqvf5bLFi8//77af78+WmjjTZaaHx+nncwS5LHr8j8tSL/wutpp52W9tprr7Tddtstdb68oeWmsbvuuqs8COXX7bnnnunf//53qjX5gJH7DNx///3p2muvLQ8se++9d/lLeY1U9xX5XOuHH35YnmNuhPpfVKUeV6SOV2YfUkvyKaDc5yKfAlzWj0+t6LZUzfJpkD/+8Y/poYceSpdddll5yviQQw4p67nR1oEbb7yx7IN3+OGHL3O+Wqj/Fv91U1LZ1yL3FVjeebEvfvGL5aMiH1T69OmTrrvuuvSzn/2spooy7ywqdthhh3LjyH+J33LLLZ8podebP/zhD2WZ5L86GqH+Wbbc5+qoo44qO7Tmg0WjbEtHH3100/9zJ9a8PFtttVXZinHAAQekRnLDDTeUrQ/L66BdC/XfYi0W66+/fmrTpk3Z3NNcfr7xxhsv8TV5/IrMXwtOPfXUdM8996SHH354uT8fv6i11lor7bTTTmnixImp1uWm3t69ey91Weqx7ityB8zRo0enE088sWHrv1KPK1LHK7MPqaVQkdeL3KF3RX8qe3nbUi3JTfu5npe2LPW6Djz66KNlx90V3SdUa/23WLBo165d2mWXXcomr4rctJufN/+rrLk8vvn8Wd7wljZ/Nct/ieRQcccdd6S//e1vqWfPniv8HrkJ8MUXXywvz6t1ue/ApEmTlros9VT3ixo+fHh5Trl///4NW/95/c8HguZ1PGPGjPLqkKXV8crsQ2olVORz5jlsfu5znwvflmpJPs2X+1gsbVnqcR2otGDm5cpXkNRF/bdkT9E///nPZa/vESNGFK+88krx7W9/u1h33XWLd955p5z+9a9/vTj33HOb5v/73/9etG3btrjyyiuLCRMmlL1h11prreLFF18sas0pp5xS9vAfM2ZMMXXq1KbHrFmzmuZZdPkvuuii4oEHHigmTZpUjBs3rjj66KOLDh06FC+//HJRa84888xy2SdPnlzW64EHHlisv/765dUx9V73zeUe7JtvvnlxzjnnLDat3up/5syZxbPPPls+8q7mqquuKv9fueLh5z//ebn933XXXcULL7xQ9ojv2bNnMXv27Kb3yD3kf/Ob33zmfUgtlcHcuXOLgQMHFptuumnx3HPPLbRf+OSTT5ZaBsvblmpl+fO0H/7wh8UTTzxRLsvo0aOLnXfeuejVq1cxZ86culgHZi5nG8imT59edOrUqbj22muX+B61WP8tGiyyXEB5x9quXbvysqEnn3yyadq+++5bXn7U3C233FL07t27nH/bbbct7r333qIW5ZVqSY98SeHSlv+0005rKquNNtqoOPTQQ4vx48cXteirX/1q0a1bt3JZNtlkk/L5xIkTG6Lum8tBIdf7a6+9tti0eqv/hx9+eInrfGUZ8yWn559/frls+UBxwAEHLFYuW2yxRRkqP+s+pJbKIB8YlrZfyK9bWhksb1uqleXPf1T169ev2GCDDco/GvJynnTSSYsFhFpeBx5ezjaQXXfddUXHjh3Ly62XpBbr38+mAwBh/FYIABBGsAAAwggWAEAYwQIACCNYAABhBAsAIIxgAQCEESwAgDCCBQAQRrAAAMIIFsAqmTZtWvlLpZdccknTuMcff7z8JcpFf6EWqH9+KwRYZaNGjUqDBw8uA8XWW2+ddtxxxzRo0KB01VVXKV1oMIIFEGLo0KFp9OjRadddd00vvvhieuaZZ1L79u2VLjQYwQIIMXv27LTddtult99+O40bNy5tv/32ShYakD4WQIhJkyalKVOmpAULFqQ333xTqUKD0mIBrLK5c+em3XffvexbkftYXH311eXpkA033FDpQoMRLIBVdtZZZ6XbbrstPf/882mdddZJ++67b+ratWu65557lC40GKdCgFUyZsyYsoXipptuSl26dEmtW7cu///oo4+ma6+9VulCg9FiAQCE0WIBAIQRLACAMIIFABBGsAAAwggWAEAYwQIACCNYAABhBAsAIIxgAQCEESwAgDCCBQAQRrAAAFKU/wcC8CzD3rZu0AAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Cell 5: Visual Diagnostics - boxplot + histogram\n", + "nBins = 25\n", + "plt.figure()\n", + "plt.hist(x, bins=nBins)\n", + "plt.title(\"Histogram\")\n", + "plt.xlabel(\"x\"); plt.ylabel(\"Counts\")\n", + "plt.show()\n", + "\n", + "plt.figure()\n", + "plt.boxplot(x, vert=False, showmeans=True)\n", + "plt.title(\"Boxplot\")\n", + "plt.xlabel(\"x\")\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv (3.14.2)", + "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.14.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/session-3-probability/probability.ipynb b/session-3-probability/probability.ipynb new file mode 100644 index 0000000..d88c41a --- /dev/null +++ b/session-3-probability/probability.ipynb @@ -0,0 +1,193 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 25, + "id": "16ca28a2", + "metadata": {}, + "outputs": [], + "source": [ + "# Week 3 — Probability Basics (Navidi Ch. 2) — Jupyter/Python Lab\n", + "# --------------------------------------------------------------\n", + "# Suggested run order: top to bottom.\n", + "# This notebook focuses on: sample spaces, events, relative frequency,\n", + "# conditional probability, independence, Bayes (via simulation), and Monte Carlo.\n", + "\n", + "import numpy as np\n", + "import math\n", + "import matplotlib.pyplot as plt\n", + "\n", + "rng = np.random.default_rng(12345)\n", + "\n", + "def relfreq(successes, trials):\n", + " return successes / trials\n", + "\n", + "def simulate_repeated(trials, sim_fn):\n", + " \"\"\"Run sim_fn() 'trials' times; sim_fn returns True/False for success.\"\"\"\n", + " x = np.fromiter((sim_fn() for _ in range(trials)), dtype=bool, count=trials)\n", + " return x.mean()\n", + "\n", + "def running_estimate(trials, sim_fn):\n", + " \"\"\"Return running estimate of P(success) over time.\"\"\"\n", + " x = np.fromiter((sim_fn() for _ in range(trials)), dtype=np.int8, count=trials)\n", + " return np.cumsum(x) / np.arange(1, trials + 1)\n", + "\n", + "def plot_running(est, p_true=None, title=\"Running estimate\"):\n", + " plt.figure()\n", + " plt.plot(est, '.k')\n", + " if p_true is not None:\n", + " plt.axhline(p_true)\n", + " plt.xlabel(\"n\")\n", + " plt.ylabel(\"estimate\")\n", + " plt.title(title)\n", + " plt.show()\n", + "\n", + "def bernoulli(p, rng=rng):\n", + " return rng.random() < p" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "id": "8e213def", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHHCAYAAABDUnkqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAARxdJREFUeJzt3Qe8FNX5//EDSBGuQBQERYoIf4kFUTomokixRGwgoJEbghQLQSkCakBsoKCigBALliQKaBSNGlARosYLKEjsBSMlKk0DCCgozP/1PfnNZnbv7N1yt85+3q/Xctnd2dnZ2dmdZ895nnMqOI7jGAAAgIComO0NAAAASCWCGwAAECgENwAAIFAIbgAAQKAQ3AAAgEAhuAEAAIFCcAMAAAKF4AYAAAQKwQ0AAAgUghvkpFNPPdVeCtXSpUtNhQoV7N+gatKkifnNb35jCtWUKVNM06ZNTaVKlUyrVq1Stt4rrrjCdOvWzeSDRx55xB7na9eujfsx33zzjalRo4Z58cUX07ptyG8EN0iJzz//3AwZMsR+WVerVs3UrFnTnHzyyeaee+4x33//PXu5QL355pvmxhtvNNu2bTO55LbbbjMLFizI2vO/9NJL5tprr7WfkYcffthuTzQKABUAuBd9tk444QRz5513mj179oQt+8UXX5gHH3zQXHfddaHbFDjocVOnTvVdv94f3b9161aTDw455BBz2WWXmd///vcpXe+XX35pLrroIlO7dm27j88991zzr3/9K67H6oeY9z1yL2eccUZKtxHxOyCBZQFfL7zwgundu7epWrWq6d+/vznuuOPM3r17zRtvvGFGjx5tPvjgA3P//fcn/OVfyE455RQbFFapUsXke3AzceJEe4LWScPrk08+MRUrZuf3lYKJXr16mfPOOy8rz//qq6/a1/7QQw/F9R7rs6WgRRQo/uUvfzGjRo0yb731lpk7d25oOf2YOPLII81pp51mgmzo0KHm3nvvtfuxS5cu5V7fzp077T7bvn27DQwrV65s7r77btO5c2ezevVqG1DFcsQRR5hJkyaF3Xb44YeXe9uQHIIblIt+Kfbt29c0btzYftEcdthhofuuvPJKs2bNGhv8JCqbJ/Xdu3eb6tWrm2zSiU8tYEGmE3ah2rx5sznwwAPjPs4POOAA8+tf/zqs66l9+/Zm3rx55q677rIn0R9//NH8+c9/tif+oPv5z39uf0SpWysVwc19991nPvvsM7NixQrTtm1be9uZZ55pn0MtZGW1rLlq1aoV9h4hu+iWQrnccccd9lePfoF6AxtXs2bNzPDhw0PXf/rpJ3PzzTebo446yp7clHehX0qRzeuROTduDsr8+fPNrbfean8l6eR/+umn2wAqMjj5+OOP42pm13PoC2zlypW2tURBjdukr+dTk32sXBE3b+Af//iHGTFihKlbt67NCTj//PPNli1bSj32V7/6lW3VateunX0N6sp77LHHYubcuNv64Ycf2l+Z2tYGDRrY9yDSunXrTM+ePe12HHrooeaaa64xixYtijuPR030v/3tb029evXs+3TssceaOXPmlFpu+vTp9j5ty89+9jPTpk0b8/jjj9v7tO/UcidqTXCb6t38imj7Ufvmd7/7nd2Pau1Rd6daAtVioZZBPY8u6tZxHCdse9T10qlTJ/tLW8FD69atzVNPPRW2jJ5j165d5tFHHw1tk3c74n3tfuI5vvV86orSNrjPr9eeaPDrfj7c/an9pmO+a9euJhWWL19uu1V00tb7q1YMHeORx5kCraOPPtrub+13teL65dCoBVeBiJbT5/eWW24x+/fvL7Xc22+/bXr06GHq1Kljl9Wxo/cjkvKK/vrXv5Y6BpKhY0RBjRvYSIsWLez3i75z4qX3X9+HyD5ablAu+nLRyVknlHior1wnFXUJjBw50n6Bqin3o48+Ms8880zMx0+ePNl+satJXk3IOrFfcskldj0u/frSyX/ChAm+wYlfgqJ+pakFSr+8dFJLxrBhw+xJV8+rL/dp06aZq666yv669lIwptc/cOBAU1xcbE+cOrnqRKwTaVn+85//2BPOBRdcYPMD9KU8ZswYc/zxx9vXIDpp6iTy9ddf28Cyfv36NuBYsmRJXK9j06ZNpkOHDvakq+1XkPG3v/3Nbu+OHTvM1VdfbZd74IEHbBCi16Ln+eGHH8y7775r34uLL77YbuOnn35qnnjiCdvEr5OVaH2x9qO2Wd1Zy5Yts12aCnLUxdWoUSP7K1rJpErIVbCngMfbLaOgTseEAiJ12ehk+/zzz5uzzz7bLvPHP/7RHocKLgcPHmxvUzCSyGsvz/Gt59dr0nHqdjXF+/mJzHMTt8tE+0fbfeKJJ/our6DfL+DX7ZHUCqvjScekjmd95hSQ6bh6/fXX7b4TdYvpefXZUcCi437WrFk28FIQ7raAbty40X4mdfIfO3asDbq1DxS8RLZode/e3e53Laf3Xet8+umnS22jtk3HlYImHQeiIPK7776La/+5x6MCLB23fgGUXqe6yLXOgw46qMz16VjX69Jxp++QQYMGmfHjx9suLmSBAyRp+/bt+snknHvuuXEtv3r1arv8ZZddFnb7qFGj7O2vvvpq6LbOnTvbi2vJkiV2mZ///OfOnj17Qrffc8899vb33nuv1LITJkyIuU16Di07e/bsUvdFW0fjxo2d4uLi0PWHH37YLtu1a1dn//79oduvueYap1KlSs62bdvCHqtlX3vttdBtmzdvdqpWreqMHDmy1GvQ38htfeyxx0K3aV/Ur1/fufDCC0O33XnnnXa5BQsWhG77/vvvnRYtWpRap5+BAwc6hx12mLN169aw2/v27evUqlXL2b17t72u9/3YY48tc11Tpkyxz/nFF1/EvR979OgRth87duzoVKhQwRk6dGjotp9++sk54ogjwo4RcbfNtXfvXue4445zunTpEnZ7jRo1wp470dde3uNbz61tiIe77JYtW+xlzZo1zm233Wb3ScuWLUPL/frXv3YOOeSQUo/Xvtfzx7po3aJ937x581Lvg177kUce6XTr1i3stkglJSWljtOrr77a3rZ8+fKw41771Ht8PPPMM/b6W2+9FXO/vPnmm3bZefPmlTqG4rm49Lp1/aabbir1HDNnzrT3ffzxx2Vuy29/+1vnxhtvdP7yl7/Y192zZ0/7uIsuuijm60B60C2FpOmXrMT6ReNySzfVdeOlX7gST27OgAEDwvIUfvnLX9q/3qoG/WpUbBJPq42o+0DrLS+1AuiXs3fb9u3bZ5vuvY455pjQdot+papZP57KjKKiorB+fe0L/br0PnbhwoW2u0otGC51f+mXZCzab0pWPeecc+z/9UvfvairQK1lq1atssvqV/W///1v++s9ldRK4t2Pyi3Rtuh2l8qn1QUWuc+8LQFq5dL2al+725yq156u4zsatcbpONFFXb3q6urYsWNYa6daINVyWNbx+fLLL5e6XHrppWHLKYFW+SdqfdM63X2gbVA3zWuvvRbqTvLub+X8aHltn44N777SvlGLmNviI3otamHzcpPO1dKm9ZXFfa3e1ii9T36v0e/icqs5/XLA3Ly3WBWf6pZXC5daK7U/n332Wft5U5eWWh+ReXRLIWkql5R4m4F1klfztr78vNQFoS+1yCDAj7ol/L7gdCJLlgKBVCQwx7ttkcu5y8bzGtT07z3xu49Vs7pL+1HdLJHLRe53P8oRUm6LugyiVbip60DUHfbKK6/YE5bWre4EnRBV3lwekftHOR/SsGHDUrdH7jOdFJXLoRN0ZJ5LKl97uo7vaHSSVRewexJWHoqOhUhl5Z80b97cNx9HuTpeCmxEXabRKNDTcaeTvrrd1GWlXCXv82sZl167gtRICuq9lNdz4YUX2i5JdTnph4oq2nRcRQYf7nN531vl/fnl/pXFDdAi8/5EXa3eZRKhoFZdt/qMKLBDZhHcoFzBjao03n///YQeF8+JJhr9YvdTnqTCRL+41BpTnm0rz2tIx+v3cn+Rq3Uo2smtZcuWoYoVlXMroFBrkVo9VHWiPAOdnJIV7TX63e593coFUWuVEsO1HTrJKd9BJ143yTlVrz1dx3c0eu2xEoWVe1OeID9yPyinKdrggmpBdPOjtH+Vi6SWJAWcev3KwfFLFo5Fj1UemVo7FMwpCV65MKpY0m3u84r7Wt3cGVGw5Q2qyqKgUw4++GAbOClHLZJ7WzIl3W4w/u233yb8WJQfwQ3KRZU/+pVbUlJiv9zKonJxfeHpl6FOjC4lceoXs+7PJfplGjn4nJIF/b4Ec4n2o5I5deL3nmgjq8r8qKtA3YwK4OKpulECZZ8+fexF+0bN8qpmGzdunG1tSMeJPhoFV3pOnRC9v/J18o3kt12JvvZcO75V3aNScJ3c3dauZLjJ1frxEms/KBBRIKjgw9vaEfm50Wt3W4S8FBz7UUuHLjqWFJiq+0rJ4UrY9g5DId59reT9eLuY3cBYrW1KyFeVViQlhKtgIt6udy+3yzRWAj3Sg5wblIvKcXWC05eOvsT9KjpUwSJnnXWW/asqIi+N0yFuNUt5JVIKHutLXvkFXgrkorXc5ArlHaiL4Lnnngs74aiJPJ4WAnULKFDwa5HzlrYrv8JLXXvKJ9JJw82X0LEhmRihWNuuoMX7/qjSxm8kYm1X5DYl8tr9ZOr4jkY/LrTvNaxBeagKSce+yur9ypq9+0H7LLLVUMMDRH5GtG/U8qIKMe96FIx5qTUmcn1u61Fkt5Fep4I4b4VhMjk3ouo25Y55AxwFXqoaU7Wdl75b1q9fH5Z7GLlteg3qHnW3CZlHyw3KRV+C+mWlX+76BeUdoVglok8++WRoDBENGa9feQoQdGJR/7q+7FQ6q371VI2qmmgpeDQK2DQgmk54GlPjn//8p20V8DaD5yKNCzNjxgzTr18/W6Kt7hmdRNzkyFitKSq3V9m4ciSUFKmARU3rShBV/oDbzK4cGzXtK8dGpa8qd9bz6iTu/tLViVKuv/5621WhbiIl7LpBTyrpeRVIqFReORrKj5k5c6bNgfHmJLnbpdfiDoCnHBa93nhfu59MHd/R/OIXv7BdU9rO8gxsp5YMlairFFyBg1pClJemgFn7Ri06bv6PWm5V2q4gQ/tKLbh6/sgRffUjSMvpvdEx6ZaCq0XH+95oX6lLUWNE6btF+XwKyvWcbvDoUoCiY6m8OTeisXr0PDqGNMyEjlMdGzqu3YRwl77n9N6640Xp2NBnTRcda+oaU6K3xgRSEvdJJ52U8PYgBdJUhYUC8+mnnzqDBg1ymjRp4lSpUsU56KCDnJNPPtmZPn2688MPP4SW+/HHH52JEyfaktLKlSs7DRs2dMaNGxe2TFml4E8++aRvmatKQJMtBY9Wzrxv3z5nzJgxTp06dZzq1avb0liV4UYrYY4sX/Ur59Zjzz77bN/t8Hu9kaXgftuqbdF6vf71r3/Z5znwwAOdunXr2jJzlalqncuWLYu5XzZt2uRceeWV9v3R+6Ry89NPP925//77Q8v84Q9/cE455RRbfqxS9qOOOsoZPXq0HSLA6+abb3YaNGjgVKxYMazsN979qPfRW6pcVjn1Qw89ZMuYtT0qfdc63cd7qbRX2679o/u82xHPa48m3uM7mVLwePzud79zmjVr5vsZUVm+n2j795133nEuuOCC0Pur90ulzYsXLw4t85///McZMGCA/YwUFRXZz4j2beR7K++++649hqtVq2aPBx0Xer+8x8SqVaucfv36OY0aNbLPeeihhzq/+tWvnLfffjtsXR999JF93CuvvOKkyoYNG5xevXo5NWvWtK9Fz/vZZ5+VWk7P6/2s6rPWu3dv+92n16bvitatW9vhJbyl9MisCvonFUESgNym7hKNVKzybf0SR/Aoz0O5Nxp4UGXbQaUEZnUZq2sqk3ldyB8EN0AAqWncWwWmnBuNXKtcCI2kiuC6/PLLbfJ4ZF5JUCjXS91ZGkMmsqsKcBHcAAGkfAmNF6NkTFXP/OlPf7LD1Cv3RvkoABBkJBQDAaQKDSWFKphRa42SPVVKq8RvAAg6Wm4AAECgMM4NAAAIFIIbAAAQKAWXc6Ph0b/66is7yBglhAAA5AeNXKOBHTXwpgabLEvBBTcKbCJnFwYAAPlhw4YN5ogjjihzmYILbtxh4bVzNKQ3AADIfZrHS40T8UxkWnDBjdsVpcCG4AYAgPwST0oJCcUAACBQCG4AAECgENwAAIBAIbgBAACBQnADAAACheAGAAAECsENAAAIFIIbAAAQKAQ3AAAgUAhuAABAoBDcAACAQCm4uaUAhPv3v/9tPvvsM1NUVGR27twZ+tu8efOYM+8CQC4iuAEKLHjxBi1Tp041Y8aMMfv37y/1uIoVK5r777/fDBw4MAtbDQDJI7gBAuyhhx4ygwcPDgte3KDl22+/Nddee23Ux+oxemzLli1p0QGQVyo4juOYArJjxw5Tq1Yts337dlOzZs1sbw6Q1habxo0b+7bKJKJChQrG+zVBiw6AXD9/k1AMBNQtt9xS7sBGIn//uC06Cp4AIBcR3AABNGXKFPOHP/whbetXgHPPPfekbf0AUB50SwEBoxaVRo0alWpxSTV1T61bt87+XwnLVFcBSCe6pYACDmxGjhyZ9sDGbb3p1auXDaS6dOliGjZsaEaPHk13FYCso1sKCFBllAKN+fPnJ/X4u+++26xYsSKhxy9fvjwskFJpuZKYtS0AkC10SwF5Pn6NuoOkrMqoSy65xJx77rmmb9++vstUqlTJrF27NjT+jXJ2yioTj7fLikEAAaQK3VJAwKllRMGMuoP0VwPtRQtsFGhMnjzZ9O7d245vo0DGS9eVfOwNRNS9pABHZeCiv0OGDAldj0XbomotAMgGWm6AgI9fM3ToUDNr1qywx69Zs8bUqFHD7Nq1yzRr1ixqC4u7rLtMoi06Wn7UqFFxLw8AqWi5YYRiIM+oBDuR8WvUuuOlICXe7qLIZdWi8/nnn8ddZq5ASF1hdE8ByCQSioE8opaUO++8M+7l1Y3UsWPHlG7DDTfcYLu64qFkY7qnAGQaLTdAjiQFu60bkbd5r6vVJt4yb3eahFS3mmh9Wq9ycPbt22cDKF2itSaplUfdWoXYPVXWpKXR7vM7JgAkyMkBM2bMcBo3buxUrVrVadeunbN8+fKoyz788MP6Zg+76HHx2r59u32M/gLZ9OCDDzoVK1a0x6P+6nrkbcXFxaHr8V4qVKjgzJ8/39mwYUNat1/rX7Jkif3r/v+SSy6Juk3Rtke3v/rqq2H3+92Wb6ZMmeL73g0ePNi5/vrrS92nfXTWWWeF3a5l83kf5Arv8ZQvx1a+bGcmtzmR83fWg5u5c+c6VapUcebMmeN88MEHzqBBg5zatWs7mzZtihrc1KxZ0/n6669Dl40bN+ZMcOO+uStWrHDmzZtnL/l0cCIzdEz4nfh0gkskkPG7jBo1KquvK9prGDJkSNjnQ3+9J3k97o477ggLCtzb8u1LX6+rvO+j9xK5DxAfHSM67qJ9xnIheIw8nt1t9m7n6aef7tx3332hz4132VR+FjYkuD4tp3Octk0/bNK9bxM5f2e9Wqp9+/ambdu2ZsaMGfa6mrY10umwYcPM2LFjSy3/yCOPmKuvvtps27atXNnWX2/5JuWzgj/66GNm2LCrSjXPq8l+xoyZpri4f0qfD/nr+uuvT8vcTOqK+uijj0yDBg1MtkybNs3m5cQzw3i8br31VjN8+PCwz5he6/TpM3Luc1XW6y+PoUMvNz3POccc1ewo+/5++eWX5vM1n4eue7n31SgqMuvWrrW3te/QPqvHRabp86XPWTxGjxptTjvtNN996SprfydD67v99jvMww/PsZ8JfTa6d+tuFr20KOZj3WVffuVl+1nQ9ZtvvtmeG8uzv37/+9+HrU/DR/i9Znfb58wpe7DOVJ/7dP4+rO4hcVVLZTW42bt3r6levbp56qmnzHnnnRe6vbi42AYvzz77rG9wc9lll9kdrTfhpJNOMrfddps59thjfZ9jz5499uLdOQqeGl4931SsWj1NrwwAAKTS/j27zYZpF8UV3GS1Wmrr1q02IbFevXpht+v6xo0bfR9z9NFHmzlz5tjA509/+pMNcDp16hR1PptJkybZlhr3osAGAAAEV1Zbbr766ivbAvPmm2+GlatqbIy///3vdt6aWH788Ufz85//3PTr1882o8XbcpPqbik102k7olWMaBTYDz/8sKCahfFf3q6UZLtl/GhdH3/8cc4eU+nqnvHz2KOPZaTbxdt070qkOyFXXDv6WjN+wviUr9fbdaMWeb99Vd7uE/c5nvvrX83s2f8bnDKbNNyB+5oiP++jRo6yXV5Lly41d0y5IyPbM3To5Wbq1Cnmpok3Zew5o0nluS+RbqmsJhTv2bPHqVSpkvPMM8+E3d6/f3+nZ8+eca+nV69eTt++fbOeUKxKF72eyGRAtxIG+ak8SXvREofLe8mXYyoymTOdl8h9Epm8HJmM6RW5rN8yySYJK9FSlWSRj1fy5XXXXRe6z5tEnYrE8lgXPWcmqsOiPXesz5U3WVV/9d5EJtrm0mXYsGF2O3Nl+4499tisb0Oqv6fyqlpKpd9XXXVV6Pq+ffucBg0aOJMmTYrr8T/99JNz9NFHO9dcc03OVEvpy0ofRJXjZqIkF5kt106EKpdS8SXhlnfruHLLr/OBtjOTX6baTzrBlHWijXwf/Zb1VmnpNUQrcY9neyJL3KN9L/iV1uuEmcsBjhuglHc7L774YhvEuJdUBzHnn3++M2vWLHtxP0e//vWvM3psJnrR9mk7tc3dunVLy48kk8LLbbfdFtq36Tr35V0puMapeeSRR5wPP/zQlo+pFNwt77700kudsWPHhpafOHGis2jRIufzzz93Vq5caVtsqlWrZsvI48E4NyhPq4uuR5b3R/sFWlZZtPfSu3fvsOXc8W3cVkD9zYdWmmgUJJT1+vX61PqayS9i7WO9P7G2rUuXLuV6HgVOmWwRyWSAk+7tysTry6XX4H4HaHv8Sv/dgNfbyqe/LVq0yOp2V8xgK3JeBTcyffp0p1GjRna8G7XkLFu2LHRf586d7Re96+qrrw4tW69ePTvo1apVq+J+LoIbxEsBS6wPdFktO2U93ntxf61H/trx/pLPdzqJuF/e+qsWrchWKO8y7r5M58lHv+bzpdvH78QW66LWFLcVOVYrRbzbmqqWmkxd4nld3n0b6yR+/PHHp3T73G5Jb2tdPJ/3yGXL836c/3+tWu53j/td5G3lcq9Htn5l+vsp74KbTCK4Qbxi/arXl51fd4bbsqMPf6wvFrVaBCF4iUc8X95+y3i7et37MpnLk2uD7nn3kV/wF+2Xf6yTd6xAIJdaOSIvCt4iWzSSeQ+8PzLc4y3yJK7rEyZMcM4444xybW+qP/ex3t9hw4Y5Q4cOLfc+yiaCmxTtHBSu8iYC67Hdu3fP6dGE81m6ErUT+bXrlwjsjsKc6X3hdxL2k0yAU95g0m1BSkdrT+QJOtOtnbH2Z2S3dLoDCm/rp/FsQzb3USoR3KRo5yC4YlVqpCoRONYXcz5+weSKaNWJmbj4tSzly3sZ64TszRPyO1mWJ3dE60tmf3u7NBVo5VJifVmtZ4kEnqkS2a20IQf2UUFOv5Bp7vQLcdXJI5AeeughM3jw4NAQ/prheuDAgaH7NSBko0aNUjYejZc7zo3GftBs2d7nReL0Xq1Zs8a88sordsDOaONMxWvChAl2fX/+85997/c7XvKRxh/SlBbRzJs3z7zxxhtm+vTpca9T+2by5Ml2zDHtQ80E7zerud6z559/3g7U2qZNG/P999+bb775JnT/P/7xD/P444+HPp/xrDNXjsMaNWqYXbt25ex2FtL5m+AGBUVfQo0bNw47CSrQWLt2bejLaMmSJaZLly4pf259US9btowvvwwGOgomFcgqGNFJ5+mnn456wtay69evt/+PPEZ0n074Gmw0KCetWAFOvDQP4AUXXJDSE7r7XhIkIOnGCafA0C1V2JTs69fsrSbjeBOJk72QX5M5ZXUVRetG8Fa6ebu88r0UvyzJjt8TrdsJSCe6pcpAt1Rhd0dp0lU/U6ZMMaNGjfJt2RHNjqtf/ZoLLdlWm3Xr1gXmV3++i6cboRBaD/Qak5lvb+jQoXbG7aDuF+QmuqVStHOQ+1/Mn332mWnevHnML9lYeTTqdpg7d67t+7/iiitK3a+uKp3kSkpKTN++fRPO7XCDJyDX6NjUfH7xcHNgRo8enfbtAspz/j6gzHuBHKHgRBOsimaBX7RoUZlJwZGPnT9/fpkJwrqvT58+Ue/XL3wFUGrB0QdsyJAhYa04Co50iRb0KHESyEUKVHSyKCv/Jh15NUA6EdwgL7qTBg0aFBaceGfXVkChQOeggw6ygY/3y9dbGVUe6rpwKYjq0aNHqW4N8WvZUcKyez+QizSrde3atc2YMWPCjl1aapCvqJZCTku0LNvbihMtfyZRkdVUsSigclt2KPlGPqGkGbmMbikEhnJqEhlvRoGMAgu1rOix5Q1sRDkGiTTFe1t2aMZHPtFxTrcTgoBuKeS0t99+O+HHqMVEgUUyj01VvgwnCQDInopZfG4gZhO5cgCSoYHcxo4dG7MLKxbyZQAg/xDcIDBdUl6xhuK/++677WjBZQU4br4MzfQAkF/olkLO0vg13qqoRJQV2Cho6dWrlw1alHzsJv8q0BkxYoS56KKLmCIBAPIYwQ1ylsaySbXI1hiSfwEgeCgFR05KVRl3ZFeU22IDAAhuKTg5N8hJ0cq4I3Nk1G2l1phYvF1RAIBgI7hBzubbRAYyClBuv/32UDCjvw888IAdYE/TKyjQSdVYNQCA/EXOTQ5P9ljI+0b5Nt5EYgU6ypVRjoymN4gcIK9jx45lrpu5nQCgcBDcZJB3nqNYkz0W8r5xW2Aiq6Q06m+0AfLKKhtnrBoAKCx0S2WwVcI7gaM7TYBuL3SR+0ZBSmSgovvUWhOrbNwPXVIAUFgIbrKYIOtOE1Do4pkDKlbri1pyRo4c6XsfXVIAUFgIbrKcIFvWCbuQ900yrS/Dhw9nHwMACG4yxR0N11vpU+hD+6s7asmSJWbu3LkxW27iaX1hHwMAhEH8snBCj6z0KfQE4lgUCKrcO979xT4GgMIexI9qqSyUOZ966qmmkPdBUVFR3IGNkoQTbeHyq6YCABQOcm4y1EqhqQS6dOli/+p6LnYPpatyS+sdPXp0aB+0b98+7mkVFNy4JeAAAMSD4KbAS8DTHXhpfY0aNTJTp04NK/WOV6wScAAAIhHcFHAJ+FtvvZXWwMsN7BIJZiJRUQYASBTBTYGWgKtFxa97KJWBVzzj1/hx9xcVZQCAZBDcpFkulieX1aKiwGLz5s0Jt9745e3EM36N3/MvW7bMrksVUkxPAQBIFMFNBugErRN1rpywo7WoKHlXAU+fPn0Syr+JlrcTz/g1kS688ELTtm1bW1FGxRMAIBmMc1OA1LqiICQy8HCDm0TGl/Fblx43duxYc+uttya8bYmOaQMAKAw7EhjnhpabAhTZVRZtFu548m+iJUzfdtttSW1briRbAwDyF8FNgVLXmOZrimytSTTx2S+vpqx1xpILydYAgPxGcJPnA+SVpwx8zJgxUYMQBSiTJk2K2T2k+2+//fZQ649ond7rrnPOOcd3HVRHAQBSieCmAEcmjlYG7qUARXkzsbZZ91977bWlgqTI67169TL33Xefb1k81VEAgFQioThNoiXaZjtZNloycTRlbXMi63LXs2jRIjtQoHJr3LL4bFePAQByHwnFOSBXRyZOdGA9bXNJSYlvYDN//vy41+W+9lwriwcABA+zgqeJm2gb2XKTzWRZBSRbtmwptV3Kj9ElWqDSt29fGzErENE67rnnHnPXXXclFCTpOd3XzqzdAIB0IuemQEYmdvN/NECf8mG8SbwjR44sM1Bx55zS5JeRk2DGS0nHjF0DAMgEcm7STC0d6o5Rq0W2Tu5+uTEKbjSCcMeOHe31eHJnkinx1mMU2IwePTrJrQcAwJBzk0sU0GR7KgG/PBtdr1u3bqiLSK1MseaBSmbsmuuuu47ABgCQUXRLFcDYN34D7el6jRo1QteVT/PEE0+k/Lk1UGCujfEDAAg2gpsCGPsmMv/Hbbnp0KFD2PN36tQp4Vm8vfwemwsVYgCAwkJwk0FqwRg8eHCoi8hN1M1Ey4ZaZhYsWBA2cnDk87tBkN/owvEENmql8Rukj+kUAACZRHBTIGPfqIWmZ8+eMSfH7NGjR1LBjUYZVtJwLlWIAQAKE+PcFMDYN26LkV9CcOTzJzrInygYevfdd03btm1tC5ECpGxXiAEAChctNwUw9o0G3fMLWBRoRT6/X/JxLAqaIru3sl0hBgAoXIxzkwWZHPtGz6WB9yJbbRTAqCtJrS1+XVju/E+JUAWYghoAAFKNuaVyXCZbNt58803f7qgRI0b4Bjbizv+kKRb8uNM1eJE4DADIFXRLBZhaYDQvVCQFJhdddFGZj1Xg1bt3b9/xcZYvX24eeOABEocBADkpJ4KbmTNnmiZNmphq1aqZ9u3bmxUrVsT1OE0foBP1eeedl/ZtzDdlJRHrtsgxbuLNEdJ1N3GY2b0BALko6zk38+bNM/379zezZ8+2gc20adPMk08+aT755BNz6KGHRn2cTqy/+MUvTNOmTc3BBx9sx3BJdZ9deYMLVR4pQdftfoq8zW+ZVFH+iwYKLIuCFe3HWM+dC/NjAQAK244Ezt9Zb7lRXsegQYPMgAEDzDHHHGODnOrVq5s5c+ZEfYwSXS+55BIzceJEG9zkwyjEkbf95je/SetIxfFUPcU7xg7VTwCAfJLV4Gbv3r1m5cqVpmvXrv/boIoV7fWSkpKoj7vppptsq466RmLZs2ePjfa8l0yPQqzrkbc9+uijaR2pOLJLSfuVJGAAQCHI6iB+W7duta0H9erVC7td1z/++GPfx7zxxhu2lWP16tVxPcekSZNsC0+2Z+COxW1FSWW3jzugngJF9T6uX7/ejB071j4XowcDAIIqr0Yo/u6778yll15qK3Xq1KkT12PGjRtny55darlp2LBhRkchdruHygpy0lVKvWjRolCrkTv/kxKCyZ8BAARVVrulFKDopL5p06aw23W9fv36pZb//PPPbQLsOeecYw444AB7eeyxx8xzzz1n/6/7I1WtWtUmHnkv6RStwijytuLi4rSXUvt1kSnYI7ABAARZVltuqlSpYlq3bm0WL14cKufWCVjXr7rqqlLLt2jRwrz33ntht91www22RUdTDKSzRSYR0eZXirztlltuSWsVUlkTdVL1BAAIqqx3S6nLSK0Ybdq0Me3atbOl4Lt27bLVU6Iy8QYNGtjcGY2Dc9xxx4U9vnbt2vZv5O3ZpuAhMoCIvM1vmWT5lZUXFRVlZaJOAAAKOrjp06eP2bJlixk/frzZuHGjadWqlVm4cGEoyVhJsIlO5FholGDtzatRF5h4u6SEJGIAQCHI+iB+mZapQfwyRS02GicnVgJzWRNlAgCQ6/JqEL9CpaBEowiXd2ybaKXnfrepuw8AgKAjuMmREYzTORKxkGsDACgUBDcZ5leeXZ7RiSNLz/2QawMAKCQENxlWVnl2stwZujVPl58nnngirqkqAAAIAoKbDPPrRkpFl5FacDRLut+6O3bsWK51AwCQTwhuMsxvBONUjE6svJ0OHTpQ+g0AKHiUgmeA3wB7uk1dUTVq1DA7d+4Muy8V5eCaAXzmzJl2qgpGIwYA5DtKwfOgMkoBh+bCUmtLeaum/PJ4NHzRFVdcUe5qLAAA8g0tN2nk16Kibigl/0q0+xJtafF7Hq9k1wsAQK6g5SYPKqNSWTWloOX222+POt5NeauxAADIJ1mfW6oQKqOiTVyZqkkt1e00ZsyY0LqUb+OdVYMB/AAAhYRqqSxVRqWqaipyUECJDGxSUY0FAEC+IOcmA9zKKLXKRAYZZd0XD81PpYRkP0yWCQAoxJwbuqUywG2pSfS+ZLu+XEyWCQAoRHRLBcCIESN8k4nJtQEAFCKCmwCMoTN16lSbZ9OjR4+Uj3wMAEC+IecmT0UbQ6ekpMTs2rUr6RweAAByETk3BSDaODkKbE499dSsbRcAANlGt1SWW19U7aS/uTK7OAAA+Y7gJsfmnEo0kZgcGwAAwpFzkyP5MomMSaNAyB24T6MR6/8Kkjp16kSeDQAgkJhbKg/zZXRdM4THasGJHJFYVVKqiurTpw8zgAMAQLdUdvjly4gCliFDhpSZg+MXGCXyeAAAgo6cmyxw55XyC3BizeAdLTCK9/EAAAQdwU2WDBw40ObYJFrxFDnhZiQqpgAAhY7gJouUPJzMzOAKjNauXWvLyKdMmULFFAAAHlRL5YBkZwbX45SDU1RUxKjEAIBA28Gs4PklmZnBveXg6tpSCxAjEwMAQMtNoOaVUlcV80kBAIKIcW4KdF4pqqQAACChOC8xrxQAANFRLZWHIsvB462yAgCgEFAtlceSrbICACDfUC2VZW6JtrqP4g06EnmMd1kqpAAACEe3VIqpRFuVTJqlW39jTYSZ6GOSWT8AAIWEbqksl2gn8hhKwAEAhWpHAoP40XKT5RLtRB5DCTgAALER3GS5RDuRx1ACDgBAbAQ3WS7RTuQxlIADABAbOTc5UqKdyGMoAQcAFJodCeTcENwAAICcxzg3AeOOa1NUVGR27tyZ0Pg5AAAUmgOyvQEom8axGTx4cFhFlRKQlaczcOBAdh8AABHolsphfuPaxDt+DgAAQcI4N3lKwcySJUvs32jj2sQ7fg4AAIWKUvAc4Tetgt+4NvGOnwMAQKEiuMkBaqnx5tXo75AhQ+z/vWPguOIZPwcAgEJFQnEOKGtaBSUN9+jRw/6/Ro0aZteuXQmNnwMAQKEhuMkBbvdT5OSZbreTAhk3mHHLwt3bAQBAOLqlckCsaRXcROOpU6eWyssBAADhKAXPIX7TKviNc+OiHBwAUCh2JDD9At1SOcTb/SRvvfVW1MDGm5dD9xQAAP9Dt1SOUotN+/btowY2Qjk4AAApDG5++ukn88orr9jckO+++87e9tVXX9m5jxI1c+ZM06RJE1OtWjV7Ql+xYkXUZZ9++mnTpk0bU7t2bVs91KpVK/PHP/7RBLE03HGcqMtQDg4AQAq7pdatW2fOOOMMs379erNnzx7TrVs3c9BBB5nbb7/dXp89e3bc65o3b54ZMWKEfYwCm2nTptnS508++cQceuihpZY/+OCDzfXXX29atGhhqlSpYp5//nkzYMAAu6weF4TAZv78+b4tNqqomjx5smnbti3l4AAApDKh+LzzzrPBjLpODjnkEPPPf/7TNG3a1CxdutQMGjQoVKocDwU0OlnPmDHDXtdJvWHDhmbYsGFm7Nixca3jpJNOMmeffba5+eabU5qQlGllJQ8rsFm2bJndVwAAFJodCZy/k+qWev31180NN9xgW0681LX05Zdfxr2evXv3mpUrV5quXbv+b4MqVrTXS0pKYj5ecdnixYttK88pp5ziu4xakrRDvJd8GKU4sgtKpeIENgAAxJZUcKMTsCp1/E7QatGJ19atW+166tWrF3a7rm/cuDHq4xS1FRUV2eBKLTbTp0+3XWN+Jk2aZCM996JWoVwUbZLMu+++287+rS4376SaAAAghcFN9+7dbW6Mq0KFCjaReMKECeass84y6aYAavXq1bZU+tZbb7U5O+oS8zNu3DgbDLmXDRs2mFzkN0mmWmx69eplFi1axOB9AACkM+dGrQdqSdBD1eKg6iX9rVOnjnnttdd8E4GjdUtVr17dPPXUUzaPx1VcXGy2bdtmnn322bjWc9lll9mgRUFAvufcaMJMtWa51VDazxqNOHJqBrXmML4NAKBQ7Ej3IH46qSqJWJVO+qtWG03weMkll5gDDzww7vWoW6l169Y2b8YNbnQS1/Wrrroq7vXoMcqtyXfeSTLdUYr9KqcYvA8AgBQHN2qd6dSpkw1mdPGOfaP7oiX3+lGXklpq1PrTrl07292lma9V3i39+/c3DRo0sLkzor9a9qijjrIBzYsvvmjHuZk1a5YJ2ijFaslR9VkkBu8DACDFwc1pp51mvv7661LdT2oq0n1+ycbR9OnTx2zZssWMHz/eJhFrUL6FCxeGkow1lo43F0WBzxVXXGG7xtRKpPFu/vSnP9n1FMJAfgzeBwBAGnJuFGxs2rTJ1K1bN+z2Tz/91Laq5Gq5da7n3EQO5Ddy5MhS9+n23r17Z2W7AAAIXM7NBRdcEKqO+s1vfmOqVq0auk+tNe+++67trkLyYs0C3rFjR3YvAACpCm4UMYkae1SO7U0eVnJwhw4dfHNEkJqB/FQ9RYUUAAApDG4efvjh0EjEo0aNshNXIjMD+Wm8GwIbAADSlHOTz3I550YtN4xpAwBAFsa5EQ28p+RWVTNpMD6vVatWJbvagqaWGc0hFTmQHy02AACkefqFe++9145Do3Ltd955x45Po9nB//Wvf5kzzzwzmVXCM5CfRh/WPFL6q+sAACDN3VIaW0bzSPXr188mFmuU4qZNm9qxar799lszY8YMk6tyuVsKAACU//ydVMuNuqLckm9VTH333Xf2/5deeql54oknklklAABASiQV3NSvX9+20EijRo3MsmXL7P+/+OKLUiPqIvnkYnVN6S8AAEhzcNOlSxfz3HPP2f8r9+aaa64x3bp1s1MgnH/++cmsEhED+alqSvtZf3UdAACkMedGY7HocsAB/y22mjt3rnnzzTdN8+bNbaWPBvTLVbmec0M5OAAAWSgF19xS3sks+/btay9Iz0B+Kgtfs2YNJeEAAMQh6XFufvjhBzuX1ObNm0udjHv27JnsagueWr8UOHr3qca7adasWcHvGwAA0hbcLFy40PTv399s3bq11H2aVFMtDUgOA/kBAJCFnBu1LnTv3t2Oa6OB/PJJrufceHNv1BWlFhtGKAYAFLodCZy/kwputFKNTHzUUUeZfJMvwQ0AAMjgIH6aoXrp0qXJPBQAACCtkmq52b17t+ndu7epW7euOf74403lypXD7v/d735nchUtNwAA5J+0l4JrioWXXnrJVKtWzbbgKInYpf/ncnADAACCLang5vrrrzcTJ040Y8eODRvvBgAAINuSikz27t1rp1ogsAEAAIEIboqLi828efNSvzUAAADZ6JbSIH133HGHWbRokWnZsmWphOK77rqrvNsFAACQueDmvffeMyeeeKL9//vvvx92nze5GAAAIC+CmyVLlqR+SwAAAFKAUicAAFCYLTcXXHCBeeSRR+zAOfp/WZ5++ulUbBsAAED6ghuNCujm0yjAIbcGAAAEZvqFfMb0CwAA5J+0T5zZpUsXs23bNt8n1n0AAADZklRwo/mkNEpxpB9++MG8/vrrqdguAACA9JeCv/vuu6H/f/jhh2bjxo1hA/stXLjQNGjQgLcCAADkR3DTqlUrm0isi1/304EHHmimT5+eyu0DAABIX3DzxRdfGOUfN23a1KxYscLUrVs3dF+VKlXMoYceaipVqpTYFgAAAGQruGncuLH9u3///lRuAwAAQHYTih999FHzwgsvhK5fe+21pnbt2qZTp05m3bp1qds6AACATAQ3t912m82vkZKSEjNjxgw7S3idOnXMNddck8wqAQAAsjdx5oYNG0yzZs3s/xcsWGB69eplBg8ebE4++WRz6qmnpmbLAAAAMtVyU1RUZL755hv7/5deesl069bN/r9atWrm+++/T2aVAAAA2Wu5UTBz2WWXmRNPPNF8+umn5qyzzrK3f/DBB6GkYwAAgLxpuZk5c6ZNHt66daudAfyQQw6xt69cudJcfPHFqd5GAACA9AY3qozq3bu3qVGjhrnxxhvNl19+aW8/6qijTOfOnZNZJQAAQPaCm7/85S/mjDPOMNWrVzfvvPOO2bNnT2jiTFVSAQAA5FVwc8stt5jZs2ebBx54wFSuXDl0u6qlVq1alcrtAwAASH9w88knn5hTTjml1O21atUy27ZtS2aVAAAA2Qtu6tevb9asWVPq9jfeeMPOOwUAAJBXwc2gQYPM8OHDzfLly+0M4V999ZX585//bEaNGmUuv/zy1G8lAABAOse5GTt2rJ088/TTTze7d++2XVRVq1a1wc2wYcOSWSUAAEBKVHAcx0n2wXv37rXdUzt37jTHHHOMHbk416miS7lB27dvNzVr1sz25gAAgBSfv5NquXFVqVLFBjUAAAB5nXMDAACQqwhuAABAoBDcAACAQMmJ4EYTcTZp0sRUq1bNtG/f3qxYsSLqshoV+Ze//KX52c9+Zi9du3Ytc3kAAFBYsh7czJs3z4wYMcJMmDDBTt1wwgknmB49epjNmzf7Lr906VLTr18/s2TJElNSUmIaNmxounfvHpq8EwAAFLZylYKnglpq2rZta2bMmGGva/wcBSwaL0fj6cSyb98+24Kjx/fv3z/m8pSCAwCQfxI5f2e15Ubj5KxcudJ2LYU2qGJFe12tMvHQIII//vijOfjgg33v14zl2iHeCwAACK6sBjdbt261LS/16tULu13XN27cGNc6xowZYw4//PCwAMlr0qRJNtJzL2oVAgAAwZX1nJvymDx5spk7d6555plnbDKyn3HjxtkmLPeyYcOGjG8nAADInHKNUFxederUMZUqVTKbNm0Ku13XNfN4WaZOnWqDm1deecW0bNky6nKa80oXAABQGLLacqPpG1q3bm0WL14cuk0JxbresWPHqI+74447zM0332wWLlxo2rRpk6GtBQAA+SCrLTeiMvDi4mIbpLRr185MmzbN7Nq1ywwYMMDerwqoBg0a2NwZuf3228348ePN448/bsfGcXNzNGlnPkzcCQAAAh7c9OnTx2zZssUGLApUWrVqZVtk3CTj9evX2woq16xZs2yVVa9evcLWo3FybrzxxoxvPwAAyC1ZH+cm0xjnBgCA/JM349wAAACkGsENAAAIFIIbAAAQKAQ3AAAgUAhuAABAoBDcAACAQCG4AQAAgUJwAwAAAoXgBgAABArBDQAACBSCGwAAECgENwAAIFAIbgAAQKAQ3AAAgEAhuAEAAIFCcAMAAAKF4AYAAAQKwQ0AAAgUghsAABAoBDcAACBQCG4AAECgENwAAIBAIbgBAACBQnADAAACheAGAAAECsENAAAIFIIbAAAQKAQ3AAAgUAhuAABAoBDcAACAQCG4AQAAgUJwAwAAAoXgBgAABArBDQAACBSCGwAAECgENwAAIFAIbgAAQKAQ3AAAgEAhuAEAAIFCcAMAAAKF4AYAAAQKwQ0AAAgUghsAABAoBDcAACBQCG4AAECgENwAAIBAIbgBAACBQnADAAACheAGAAAECsENAAAIFIIbAAAQKAQ3AAAgUAhuAABAoGQ9uJk5c6Zp0qSJqVatmmnfvr1ZsWJF1GU/+OADc+GFF9rlK1SoYKZNm5bRbQUAALkvq8HNvHnzzIgRI8yECRPMqlWrzAknnGB69OhhNm/e7Lv87t27TdOmTc3kyZNN/fr1M769AAAg92U1uLnrrrvMoEGDzIABA8wxxxxjZs+ebapXr27mzJnju3zbtm3NlClTTN++fU3VqlUzvr0AACD3ZS242bt3r1m5cqXp2rXr/zamYkV7vaSkJGXPs2fPHrNjx46wCwAACK6sBTdbt241+/btM/Xq1Qu7Xdc3btyYsueZNGmSqVWrVujSsGHDlK0bAADknqwnFKfbuHHjzPbt20OXDRs2ZHuTAABAGh1gsqROnTqmUqVKZtOmTWG363oqk4WVm0N+DgAAhSNrLTdVqlQxrVu3NosXLw7dtn//fnu9Y8eO2dosAACQ57LWciMqAy8uLjZt2rQx7dq1s+PW7Nq1y1ZPSf/+/U2DBg1s3oybhPzhhx+G/v/ll1+a1atXm6KiItOsWbNsvhQAAJAjshrc9OnTx2zZssWMHz/eJhG3atXKLFy4MJRkvH79eltB5frqq6/MiSeeGLo+depUe+ncubNZunRpVl4DAADILRUcx3FMAVEpuKqmlFxcs2bNbG8OAABI8fk78NVSAACgsBDcAACAQCG4AQAAgUJwAwAAAoXgBgAABArBDQAACBSCGwAAECgENwAAIFAIbgAAQKAQ3AAAgEAhuAEAAIFCcAMAAAKF4AYAAAQKwQ0AAAgUghsAABAoBDcAACBQCG4AAECgENwAAIBAIbgBAACBQnADAAACheAGAAAECsENAAAIFIIbAAAQKAQ3AAAgUAhuAABAoBDcAACAQCG4AQAAgUJwAwAAAoXgBgAABArBDQAACBSCGwAAECgENwAAIFAIbgAAQKAQ3AAAgEAhuAEAAIFCcAMAAAKF4AYAAAQKwQ0AAAgUghsAABAoBDcAACBQCG4AAECgENwAAIBAIbgBAACBQnADAAACheAGAAAECsENAAAIFIIbAAAQKAQ3AAAgUAhuAABAoBDcAACAQCG4AQAAgUJwAwAAAoXgBgAABEpOBDczZ840TZo0MdWqVTPt27c3K1asKHP5J5980rRo0cIuf/zxx5sXX3wxY9sKAAByW9aDm3nz5pkRI0aYCRMmmFWrVpkTTjjB9OjRw2zevNl3+TfffNP069fPDBw40LzzzjvmvPPOs5f333/fZNu///1vs2TJEvsXAABkRwXHcRyTRWqpadu2rZkxY4a9vn//ftOwYUMzbNgwM3bs2FLL9+nTx+zatcs8//zzods6dOhgWrVqZWbPnh3z+Xbs2GFq1apltm/fbmrWrJmy1/HQQw+ZwYMH2+2vWLGiuf/++20ABgAAyi+R83dWW2727t1rVq5cabp27fq/DapY0V4vKSnxfYxu9y4vaumJtvyePXvsDvFeUk0tNW5gI/o7ZMgQWnAAAMiCrAY3W7duNfv27TP16tULu13XN27c6PsY3Z7I8pMmTbKRnntRq1CqffbZZ6HAxqXXtWbNmpQ/FwAAyPGcm3QbN26cbcJyLxs2bEj5czRv3ty2OHlVqlTJNGvWLOXPBQAAcji4qVOnjg0CNm3aFHa7rtevX9/3Mbo9keWrVq1q++a8l1Q74ogjbI6NXovo7x/+8Ad7OwAAKKDgpkqVKqZ169Zm8eLFodvUvaPrHTt29H2MbvcuLy+//HLU5TNFycNr16611VL6SzIxAADZcYDJMpWBFxcXmzZt2ph27dqZadOm2WqoAQMG2Pv79+9vGjRoYHNnZPjw4aZz587mzjvvNGeffbaZO3euefvtt23LSbappYbWGgAACjy4UWn3li1bzPjx421SsEq6Fy5cGEoaXr9+fVg+S6dOnczjjz9ubrjhBnPdddfZfJcFCxaY4447LouvAgAA5Iqsj3OTaeka5wYAAKRP3oxzAwAAkGoENwAAIFAIbgAAQKAQ3AAAgEAhuAEAAIFCcAMAAAKF4AYAAAQKwQ0AAAgUghsAABAoWZ9+IdPcAZk10iEAAMgP7nk7nokVCi64+e677+zfhg0bZntTAABAEudxTcNQloKbW2r//v3mq6++MgcddJCpUKFCyqNKBU0bNmxg3qo0Yj9nDvua/RwkHM/5va8VriiwOfzww8Mm1PZTcC032iFHHHFEWp9DbySTcqYf+zlz2Nfs5yDheM7ffR2rxcZFQjEAAAgUghsAABAoBDcpVLVqVTNhwgT7F+nDfs4c9jX7OUg4ngtnXxdcQjEAAAg2Wm4AAECgENwAAIBAIbgBAACBQnADAAACheAmRWbOnGmaNGliqlWrZtq3b29WrFiRqlUXjNdee82cc845dvRJjR69YMGCsPuV+z5+/Hhz2GGHmQMPPNB07drVfPbZZ2HLfPvtt+aSSy6xg0bVrl3bDBw40OzcuTPDryR3TZo0ybRt29aO0H3ooYea8847z3zyySdhy/zwww/myiuvNIcccogpKioyF154odm0aVPYMuvXrzdnn322qV69ul3P6NGjzU8//ZThV5PbZs2aZVq2bBkaxKxjx47mb3/7W+h+9nN6TJ482X5/XH311ezrFLvxxhvtvvVeWrRokZvHtKqlUD5z5851qlSp4syZM8f54IMPnEGDBjm1a9d2Nm3axK5NwIsvvuhcf/31ztNPP60KPueZZ54Ju3/y5MlOrVq1nAULFjj//Oc/nZ49ezpHHnmk8/3334eWOeOMM5wTTjjBWbZsmfP66687zZo1c/r168f78H969OjhPPzww87777/vrF692jnrrLOcRo0aOTt37gzto6FDhzoNGzZ0Fi9e7Lz99ttOhw4dnE6dOoXu/+mnn5zjjjvO6dq1q/POO+/Y961OnTrOuHHj2M8ezz33nPPCCy84n376qfPJJ5841113nVO5cmW779nP6bFixQqnSZMmTsuWLZ3hw4dzTKfYhAkTnGOPPdb5+uuvQ5ctW7bk5HcHwU0KtGvXzrnyyitD1/ft2+ccfvjhzqRJk1Kx+oIUGdzs37/fqV+/vjNlypTQbdu2bXOqVq3qPPHEE/b6hx9+aB/31ltvhZb529/+5lSoUMH58ssvM/wK8sPmzZvtPvv73/8e2qc6AT/55JOhZT766CO7TElJib2uL6SKFSs6GzduDC0za9Ysp2bNms6ePXuy8Cryx89+9jPnwQcfZD+nwXfffec0b97cefnll53OnTuHghuO6dQGN/rx6CfX9jPdUuW0d+9es3LlSttF4p2/StdLSkrKu3r8ny+++MJs3LgxbD9rjhF1Abr7WX/VFdWmTZvQMlpe78fy5cvZlz62b99u/x588MH2r47lH3/8MWw/q9m5UaNGYfv5+OOPN/Xq1Qst06NHDztR3gcffMB+9rFv3z4zd+5cs2vXLts9xX5OPXWHqLvDe+xyTKeeUgGUOtC0aVObAqBuJsm1Y7rgJs5Mta1bt9ovLu+bJbr+8ccfZ227gkaBjfjtZ/c+/VUfrtcBBxxgT9zuMvif/fv327yEk08+2Rx33HGhfVilShUbJJa1n/3eB+/7hP967733bDCjXATlIDzzzDPmmGOOMatXr2Y/p5ACx1WrVpm33nqr1H0c06mjH5OPPPKIOfroo83XX39tJk6caH75y1+a999/P+f2M8ENUMC/dPWl9MYbb2R7UwJLJwEFMmohe+qpp0xxcbH5+9//nu3NCpQNGzaY4cOHm5dfftkWdCB9zjzzzND/lSyvYKdx48Zm/vz5tsgjl9AtVU516tQxlSpVKpURruv169cv7+rxf9x9WdZ+1t/NmzeH3a8sfFVQ8V6Eu+qqq8zzzz9vlixZYo444ojQ7dpP6mrdtm1bmfvZ733wvk/4L/2SbdasmWndurWtVDvhhBPMPffcw35OIXWH6HN/0kkn2ZZaXRRA3nvvvfb/ahngmE4PtdL8v//3/8yaNWty7pgmuEnBl5e+uBYvXhzW3K/rao5Gahx55JH24PfuZ/XTKpfG3c/6qw+Wvuxcr776qn0/9AsD/y2nV2Cj7hHtG+1XLx3LlStXDtvPKhVXv7p3P6u7xRtI6lezyp3V5YLodCzu2bOH/ZxCp59+uj0e1ULmXpR3p3wQ9/8c0+mhYTY+//xzOzxHzn13pDQ9uYBLwVW188gjj9iKncGDB9tScG9GOOKrdlB5oC46NO+66y77/3Xr1oVKwbVfn332Wefdd991zj33XN9S8BNPPNFZvny588Ybb9jqCUrB/+fyyy+35fRLly4NK+fcvXt3WDmnysNfffVVW87ZsWNHe4ks5+zevbstJ1+4cKFTt25dSsEjjB071lahffHFF/Z41XVV7r300kvs5zTzVktxTKfOyJEj7XeHjul//OMftqRbpdyqusy17w6CmxSZPn26fVM13o1KwzXOChKzZMkSG9REXoqLi0Pl4L///e+devXq2WDy9NNPt+OHeH3zzTc2mCkqKrLlhQMGDLBBE/7Lb//qorFvXAoWr7jiClu2XL16def888+3AZDX2rVrnTPPPNM58MAD7ZebvvR+/PFHdrPHb3/7W6dx48b2O0Ff4Dpe3cCG/ZzZ4IZjOjX69OnjHHbYYfaYbtCggb2+Zs2anNzPFfRPatuCAAAAsoecGwAAECgENwAAIFAIbgAAQKAQ3AAAgEAhuAEAAIFCcAMAAAKF4AYAAAQKwQ0AAAgUghsAABAoBDcAACBQDsj2BgBAeZ166qmmZcuWplq1aubBBx80VapUMUOHDjU33ngjOxcoQLTcAAiERx991NSoUcMsX77c3HHHHeamm24yL7/8crY3C0AWMHEmgEC03Ozbt8+8/vrrodvatWtnunTpYiZPnpzVbQOQebTcAAgEdUt5HXbYYWbz5s1Z2x4A2UNwAyAQKleuHHa9QoUKZv/+/VnbHgDZQ3ADAAACheAGAAAECsENAAAIFKqlAABAoNByAwAAAoXgBgAABArBDQAACBSCGwAAECgENwAAIFAIbgAAQKAQ3AAAgEAhuAEAAIFCcAMAAAKF4AYAAAQKwQ0AAAgUghsAAGCC5P8DNG9CrnNpi9sAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Final estimate: 0.5\n" + ] + } + ], + "source": [ + "# Cell 1 — Basic: coin flips as Bernoulli trials\n", + "p = 0.5\n", + "N = 500\n", + "\n", + "est = running_estimate(N, lambda: bernoulli(p))\n", + "plot_running(est, p_true=p, title=\"Coin: running estimate of P(Heads)=0.5\")\n", + "\n", + "print(\"Final estimate:\", est[-1])" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "c7e376bc", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Simulated P(Red): 0.35156\n", + "Exact P(Red): 0.35\n" + ] + } + ], + "source": [ + "# Cell 2 — Law of Total Probability via simulation (urn example)\n", + "# Urn A chosen with probability 0.3; Urn B with probability 0.7\n", + "# Urn A: 70% red; Urn B: 20% red\n", + "pA = 0.3\n", + "p_red_A = 0.7\n", + "p_red_B = 0.2\n", + "\n", + "N = 200_000\n", + "\n", + "def draw_red():\n", + " # choose urn then draw color\n", + " choose_A = bernoulli(pA)\n", + " if choose_A:\n", + " return bernoulli(p_red_A)\n", + " else:\n", + " return bernoulli(p_red_B)\n", + "\n", + "p_hat = simulate_repeated(N, draw_red)\n", + "\n", + "p_exact = pA * p_red_A + (1 - pA) * p_red_B\n", + "print(\"Simulated P(Red):\", p_hat)\n", + "print(\"Exact P(Red): \", p_exact)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "cca32796", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0, 0, 0, 1, 0, 1, 0, 0, 1, 1], dtype=int8)" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "trials = 10\n", + "sim_fn = lambda: rng.random() < 0.5\n", + "x = np.fromiter((sim_fn() for _ in range(trials)), dtype=np.int8, count=trials)\n", + "# x might be: array([1, 0, 1, 1, 0], dtype=int8)\n", + "x" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "d3782b57", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([1, 0, 1, 1, 0, 1, 1, 1, 1, 1], dtype=int8)" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x = np.fromiter((sim_fn() for _ in range(trials)), dtype=np.int8, count=trials)\n", + "x" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv (3.14.2)", + "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.14.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}