{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "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": 2,
   "id": "8e213def",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHHCAYAAABDUnkqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAR7hJREFUeJzt3Qe4FNX9//FzQTqCkS5dIXZQERBMLIiiGAs2VBSCSEnEhoggCWCJEImIQRRUFEsiYEMTDUoAC6EpaIwdIgIiVQUEBATm/3zO7z+buXtn926ZbXPfr+cZlp2dnXJ27s53z/meM0WO4zgGAAAgJMrlegcAAACCRHADAABCheAGAACECsENAAAIFYIbAAAQKgQ3AAAgVAhuAABAqBDcAACAUCG4AQAAoUJwg7x02mmn2amsevPNN01RUZF9DKtmzZqZX//616asGjt2rDn00ENN+fLlzXHHHRfYen/729+aM8880xSCqVOn2vP8q6++Svg93377ralWrZp57bXXMrpvKGwENwjEf//7X9O/f3/7ZV25cmVTo0YNc/LJJ5sHHnjA/Pjjj5RyGbVgwQIzatQos2XLFpNP7rnnHjNz5sycbf+NN94wQ4YMsX8jTzzxhN2fWBQAKgBwJ/1ttW7d2tx3331m9+7dxZZduXKleeyxx8ztt98emafAQe/705/+5Lt+fT56ffPmzaYQ1KpVy1x77bXm97//faDrXbt2rbnsssvMQQcdZMv4ggsuMF9++WVC79UPMe9n5E5nn312oPuIxB2QxLKAr1dffdVceumlplKlSqZnz57mmGOOMXv27DHz5883t956q/n444/NI488kvSXf1l2yimn2KCwYsWKptCDmzvuuMNeoHXR8Pr8889NuXK5+X2lYOKSSy4xF154YU62P3fuXHvsU6ZMSegz1t+WghZRoPjCCy+YwYMHm3fffddMmzYtspx+TDRv3tycfvrpJswGDBhg/vznP9ty7NSpU9rr2759uy2zrVu32sCwQoUK5v777zennnqq+eCDD2xAVZpGjRqZ0aNHF5t3yCGHpL1vSA3BDdKiX4qXX365adq0qf2iadCgQeS16667zqxYscIGP8nK5UV9586dpmrVqiaXdOFTDViY6YJdVm3cuNFUqVIl4fP8gAMOMFdddVWxpqf27dub6dOnm3HjxtmL6E8//WT+8pe/2At/2B155JH2R5SatYIIbh566CGzfPlys2TJEtO2bVs775xzzrHbUA1ZvJo1V82aNYt9RsgtmqWQlnvvvdf+6tEvUG9g42rRooW58cYbI8/37t1r7rrrLnPYYYfZi5vyLvRLKbp6PTrnxs1BmTFjhvnDH/5gfyXp4n/GGWfYACo6OPnss88SqmbXNvQFtnTpUltboqDGrdLX9lRlX1quiJs38K9//csMGjTI1KlTx+YEdOvWzWzatKnEe3/1q1/ZWq127drZY1BT3lNPPVVqzo27r5988on9lal9bdiwof0Moq1atcqcf/75dj/q1q1rbr75ZvP6668nnMejKvprrrnG1KtXz35ORx99tHn88cdLLDdhwgT7mvblZz/7mTnxxBPNX//6V/uayk41d6LaBLeq3s2viFWOKpsbbrjBlqNqe9TcqZpA1VioZlDb0aRmHcdxiu2Pml46duxof2kreGjTpo15/vnniy2jbezYscM8+eSTkX3y7keix+4nkfNb21NTlPbB3b6OPdng1/37cMtT5aZzvnPnziYIixcvts0qumjr81Uths7x6PNMgdbhhx9uy1vlrlpcvxwa1eAqENFy+vu9++67zf79+0ss995775kuXbqY2rVr22V17ujziKa8or/97W8lzoFU6BxRUOMGNnLEEUfY7xd95yRKn7++D5F71NwgLfpy0cVZF5REqK1cFxU1Cdxyyy32C1RVuZ9++ql56aWXSn3/mDFj7Be7quRVhawLe48ePex6XPr1pYv/yJEjfYMTvwRF/UpTDZR+eemilorrr7/eXnS1XX25jx8/3gwcOND+uvZSMKbj79Onj+nVq5e9cOriqguxLqTxfP/99/aCc9FFF9n8AH0p33bbbebYY4+1xyC6aOoism7dOhtY1q9f3wYc8+bNS+g4NmzYYE466SR70dX+K8j4xz/+Yfd327Zt5qabbrLLPfroozYI0bFoO7t27TIffvih/SyuvPJKu49ffPGFefbZZ20Vvy5WovWVVo7aZzVnLVq0yDZpKshRE1eTJk3sr2glkyohV8GeAh5vs4yCOp0TCojUZKOL7d///ndz7rnn2mWefvppex4quOzXr5+dp2AkmWNP5/zW9nVMOk/dpqZE/36i89zEbTJR+Wi/jz/+eN/lFfT7BfyaH021sDqfdE7qfNbfnAIynVfvvPOOLTtRs5i2q78dBSw67x9++GEbeCkId2tA169fb/8mdfEfOnSoDbpVBgpeomu0zjrrLFvuWk6fu9b54osvlthH7ZvOKwVNOg9EQeQPP/yQUPm556MCLJ23fgGUjlNN5FrngQceGHd9Otd1XDrv9B3St29fM2LECNvEhRxwgBRt3bpVP5mcCy64IKHlP/jgA7v8tddeW2z+4MGD7fy5c+dG5p166ql2cs2bN88uc+SRRzq7d++OzH/ggQfs/P/85z8llh05cmSp+6RtaNlJkyaVeC3WOpo2ber06tUr8vyJJ56wy3bu3NnZv39/ZP7NN9/slC9f3tmyZUux92rZt99+OzJv48aNTqVKlZxbbrmlxDHoMXpfn3rqqcg8lUX9+vWdiy++ODLvvvvus8vNnDkzMu/HH390jjjiiBLr9NOnTx+nQYMGzubNm4vNv/zyy52aNWs6O3futM/1uR999NFx1zV27Fi7zZUrVyZcjl26dClWjh06dHCKioqcAQMGRObt3bvXadSoUbFzRNx9c+3Zs8c55phjnE6dOhWbX61atWLbTvbY0z2/tW3tQyLcZTdt2mSnFStWOPfcc48tk1atWkWWu+qqq5xatWqVeL/KXtsvbdK6RWXfsmXLEp+Djr158+bOmWeeWWxetIULF5Y4T2+66SY7b/HixcXOe5Wp9/x46aWX7PN333231HJZsGCBXXb69OklzqFEJpeOW8/vvPPOEtuYOHGife2zzz6Luy/XXHONM2rUKOeFF16wx33++efb91122WWlHgcyg2YppEy/ZKW0XzQut+ummm689AtXEsnN6d27d7E8hV/+8pf20durQb8aFZskUmsjaj7QetOlWgD9cvbu2759+2zVvddRRx0V2W/Rr1RV6yfSM6N69erF2vVVFvp16X3vrFmzbHOVajBcav7SL8nSqNyUrHreeefZ/+uXvjupqUC1ZcuWLbPL6lf1119/bX+9B0m1JN5yVG6J9kXzXeo+rSaw6DLz1gSolkv7q7J29zmoY8/U+R2LauN0nmhSU6+aujp06FCstlM1kKo5jHd+zp49u8R09dVXF1tOCbTKP1Htm9bploH2Qc00b7/9dqQ5yVveyvnR8to/nRveslLZqEbMrfERHYtq2LzcpHPVtGl98bjH6q2N0ufkd4x+k8vtzemXA+bmvZXW41PN8qrhUm2lyvPll1+2f29q0lLtI7KPZimkTN0lJdFqYF3kVb2tLz8vNUHoSy06CPCjZgm/LzhdyFKlQCCIBOZE9y16OXfZRI5BVf/eC7/7XlWru1SOamaJXi663P0oR0i5LWoyiNXDTU0Houawf/7zn/aCpXWrOUEXRHVvTkd0+SjnQxo3blxifnSZ6aKoXA5doKPzXII89kyd37HoIqsmYPcirDwUnQvR4uWftGzZ0jcfR7k6XgpsRE2msSjQ03mni76a3dRkpVwl7/a1jEvHriA1moJ6L+X1XHzxxbZJUk1O+qGiHm06r6KDD3db3s9WeX9+uX/xuAFadN6fqKnVu0wyFNSq6VZ/IwrskF0EN0gruFEvjY8++iip9yVyoYlFv9j9pJNUmOwXl2pj0tm3dI4hE8fv5f4iV+1QrItbq1atIj1W1J1bAYVqi1TroV4nyjPQxSlVsY7Rb773uJULotoqJYZrP3SRU76DLrxuknNQx56p8zsWHXtpicLKvUknyI8uB+U0xRpcUDWIbn6Uyle5SKpJUsCp41cOjl+ycGn0XuWRqbZDwZyS4JULox5LmuduV9xjdXNnRMGWN6iKR0GnHHzwwTZwUo5aNHdeKl263WD8u+++S/q9SB/BDdKinj/6lbtw4UL75RaPuovrC0+/DHVhdCmJU7+Y9Xo+0S/T6MHnlCzo9yWYT1SOSubUhd97oY3uVeZHTQVqZlQAl0ivGyVQdu/e3U4qG1XLqzfbsGHDbG1DJi70sSi40jZ1QfT+ytfFN5rffiV77Pl2fqt3j7qC6+Lu1nalwk2u1o+X0spBgYgCQQUf3tqO6L8bHbtbI+Sl4NiPajo06VxSYKrmKyWHK2HbOwyFeMtayfuJNjG7gbFq25SQr15a0ZQQrg4TiTa9e7lNpqUl0CMzyLlBWtQdVxc4fenoS9yvR4d6sEjXrl3to3oReWmcDnF7s6Qrma7gpX3JK7/AS4FcrJqbfKG8AzURvPLKK8UuOKoiT6SGQM0CChT8auS8XduVX+Glpj3lE+mi4eZL6NyQbIxQrH1X0OL9fNTTxm8kYu1X9D4lc+x+snV+x6IfFyp7DWuQDvVC0rmvbvV+3Zq95aAyi6411PAA0X8jKhvVvKiHmHc9Csa8VBsTvT639ii62UjHqSDO28MwlZwbUe825Y55AxwFXuo1pt52XvpuWb16dbHcw+h90zGoedTdJ2QfNTdIi74E9ctKv9z1C8o7QrG6iD733HORMUQ0ZLx+5SlA0IVF7ev6slPXWbWrBzWqarJdwWNRwKYB0XTB05ga//73v22tgLcaPB9pXJgHH3zQXHHFFbaLtppndBFxkyNLq01Rd3t1G1eOhJIiFbCoal0JosofcKvZlWOjqn3l2Kjrq7o7a7u6iLu/dHWhlOHDh9umCjUTKWHXDXqCpO0qkFBXeeVoKD9m4sSJNgfGm5Pk7peOxR0ATzksOt5Ej91Pts7vWH7xi1/YpintZzoD26kmQ13U1RVcgYNqQpSXpoBZZaMaHTf/RzW36tquIENlpRpcbT96RF/9CNJy+mx0TrpdwVWj4/1sVFZqUtQYUfpuUT6fgnJt0w0eXQpQdC6lm3MjGqtH29E5pGEmdJ7q3NB57SaEu/Q9p8/WHS9K54b+1jTpXFPTmBK9NSaQkrhPOOGEpPcHAchQLyyUMV988YXTt29fp1mzZk7FihWdAw880Dn55JOdCRMmOLt27Yos99NPPzl33HGH7VJaoUIFp3Hjxs6wYcOKLROvK/hzzz3n281VXUBT7Qoeqzvzvn37nNtuu82pXbu2U7VqVds1Vt1wY3Vhju6+6tedW+8999xzfffD73iju4L77av2Rev1+vLLL+12qlSp4tSpU8d2M1c3Va1z0aJFpZbLhg0bnOuuu85+Pvqc1N38jDPOcB555JHIMpMnT3ZOOeUU2/1YXdkPO+ww59Zbb7VDBHjdddddTsOGDZ1y5coV6/abaDnqc/R2VY7XnXrKlCm2G7P2R13ftU73/V7q2qt9V/noNe9+JHLssSR6fqfSFTwRN9xwg9OiRQvfvxF1y/cTq3zff/9956KLLop8vvq81LV5zpw5kWW+//57p3fv3vZvpHr16vZvRGUb/dnKhx9+aM/hypUr2/NB54U+L+85sWzZMueKK65wmjRpYrdZt25d51e/+pXz3nvvFVvXp59+at/3z3/+0wnKmjVrnEsuucSpUaOGPRZtd/ny5SWW03a9f6v6W7v00kvtd5+OTd8Vbdq0scNLeLvSI7uK9E8QQRKA/KbmEo1UrO7b+iWO8FGeh3JvNPCgum2HlRKY1WSspqls5nWhcBDcACGkqnFvLzDl3GjkWuVCaCRVhNdvfvMbmzwenVcSFsr1UnOWxpCJbqoCXAQ3QAgpX0LjxSgZU71nnnnmGTtMvXJvlI8CAGFGQjEQQuqhoaRQBTOqrVGyp7rSKvEbAMKOmhsAABAqjHMDAABCheAGAACESpnLudHw6N98840dZIwuhAAAFAaNXKOBHTXwpgabjKfMBTcKbKLvLgwAAArDmjVrTKNGjeIuU+aCG3dYeBWOhvQGAAD5T/fxUuVEIjcyLXPBjdsUpcCG4AYAgMKSSEoJCcUAACBUCG4AAECoENwAAIBQIbgBAAChQnADAABCheAGAACECsENAAAIFYIbAAAQKgQ3AAAgVAhuAABAqOQ0uHn77bfNeeedZ+/wqeGUZ86cWep73nzzTXPCCSeYSpUqmRYtWpipU6dmZV8BAEBhyGlws2PHDtO6dWszceLEhJZfuXKlOffcc83pp59uPvjgA3PTTTeZa6+91rz++usmH3z99ddmxowZdtL/AQBA9uX0xpnnnHOOnRI1adIk07x5c3PffffZ50ceeaSZP3++uf/++02XLl1MLk2ZMsX07dvXOI5jn6sm6tFHHzV9+vTJ6X4BAFDWFFTOzcKFC03nzp2LzVNQo/mx7N69294m3TsFTbU03sBG9P/+/ftTgwMAQJYVVHCzfv16U69evWLz9FwBy48//uj7ntGjR5uaNWtGpsaNGwe+X8uXLy8W2Lj27dtnVqxYEfj2AABASIKbVAwbNsxs3bo1Mq1ZsybwbbRs2dI2Q0UrX768TXoGAADZU1DBTf369c2GDRuKzdPzGjVqmCpVqvi+R72q9Lp3ClqjRo1sfo03wClXrpyZPHmyfQ0AAJSRhOJkdejQwbz22mvF5s2ePdvOzzUlDnvzf7RPBDYAAJSx4Gb79u3FclLU1VtdvA8++GDTpEkT26S0du1a89RTT9nXBwwYYB588EEzZMgQc80115i5c+fabtevvvqqyQcKZi699NJc7wYAAGVaTpul3nvvPXP88cfbSQYNGmT/P2LECPt83bp1ZvXq1ZHl1Q1cgYxqazQ+jrqEP/bYYznvBg4AAPJHkePXzSfE1LNKvaaUXJyJ/BsAAJDb63dBJRQDAACUhuAGAACECsENAAAIFYIbAAAQKgQ3AAAgVAhuAABAqBDcAACAUCG4AQAAoUJwAwAAQoXgBgAAhArBDQAACBWCGwAAECoENwAAIFQIbgAAQKgQ3AAAgFAhuAEAAKFCcAMAAEKF4AYAAIQKwQ0AAAgVghsAABAqBDcAACBUCG4AAECoENwAAIBQIbgBAAChQnADAABCheAGAACECsENAAAIFYIbAAAQKgQ3AAAgVAhuAABAqBDcAACAUCG4AQAAoUJwAwAAQoXgBgAAhArBDQAACBWCGwAAECoENwAAIFQIbgAAQKgQ3AAAgFAhuAEAAKFCcAMAAEKF4AYAAIQKwQ0AAAgVghsAABAqBDcAACBUCG4AAECoENwAAIBQIbgBAAChQnADAABCheAGAACECsENAAAIlZwHNxMnTjTNmjUzlStXNu3btzdLliyJuexPP/1k7rzzTnPYYYfZ5Vu3bm1mzZqV1f0FAAD5LafBzfTp082gQYPMyJEjzbJly2yw0qVLF7Nx40bf5X/3u9+ZyZMnmwkTJphPPvnEDBgwwHTr1s28//77Wd93AACQn4ocx3FytXHV1LRt29Y8+OCD9vn+/ftN48aNzfXXX2+GDh1aYvlDDjnEDB8+3Fx33XWReRdffLGpUqWKeeaZZxLa5rZt20zNmjXN1q1bTY0aNQI8GgAAkCnJXL9zVnOzZ88es3TpUtO5c+f/7Uy5cvb5woULfd+ze/du2xzlpcBm/vz5Mbej96hAvBMAAAivnAU3mzdvNvv27TP16tUrNl/P169f7/seNVmNGzfOLF++3NbyzJ4927z44otm3bp1MbczevRoG+m5k2qGAABAeOU8oTgZDzzwgGnZsqU54ogjTMWKFc3AgQNN7969bY1PLMOGDbNVWO60Zs2arO4zAAAoI8FN7dq1Tfny5c2GDRuKzdfz+vXr+76nTp06ZubMmWbHjh1m1apV5rPPPjPVq1c3hx56aMztVKpUybbNeScAABBeOQtuVPPSpk0bM2fOnMg8NTXpeYcOHeK+V3k3DRs2NHv37jUvvPCCueCCC7KwxwAAoBAckMuNqxt4r169zIknnmjatWtnxo8fb2tl1NQkPXv2tEGM8mZk8eLFZu3atea4446zj6NGjbIB0ZAhQ3J5GAAAII/kNLjp3r272bRpkxkxYoRNIlbQokH53CTj1atXF8un2bVrlx3r5ssvv7TNUV27djVPP/20Oeigg3J4FAAAIJ/kdJybXGCcGwAACk9BjHMDAACQCQQ3AAAgVAhuAABAqBDcAACAUCG4AQAAoUJwAwAAQoXgBgAAhArBDQAACBWCGwAAECoENwAAIFQIbgAAQKgQ3AAAgFAhuAEAAKFCcAMAAEKF4AYAAIQKwQ0AAAgVghsAABAqBDcAACBUCG4AAECoENwAAIBQIbgBAAChQnADAABCheAGAACECsENAAAIFYIbAAAQKgQ3AAAgVAhuAABAqBDcAACAUCG4AQAAoUJwAwAAQoXgBgAAhArBDQAACBWCGwAAECoENwAAIFQIbgAAQKgQ3AAAgFAhuAEAAKFCcAMAAEKF4AYAAIQKwQ0AAAgVghsAABAqBDcAACBUCG4AAECoENwAAIBQIbgBAAChQnADAABCheAGAACECsENAAAIFYIbAAAQKgQ3AAAgVAhuAABAqOQ8uJk4caJp1qyZqVy5smnfvr1ZsmRJ3OXHjx9vDj/8cFOlShXTuHFjc/PNN5tdu3ZlbX8BAEB+y2lwM336dDNo0CAzcuRIs2zZMtO6dWvTpUsXs3HjRt/l//rXv5qhQ4fa5T/99FMzZcoUu47bb7896/sOAADyU06Dm3Hjxpm+ffua3r17m6OOOspMmjTJVK1a1Tz++OO+yy9YsMCcfPLJ5sorr7S1PWeddZa54oorSq3tAQAAZUfOgps9e/aYpUuXms6dO/9vZ8qVs88XLlzo+56OHTva97jBzJdffmlee+0107Vr16ztNwAAyG8H5GrDmzdvNvv27TP16tUrNl/PP/vsM9/3qMZG7/vFL35hHMcxe/fuNQMGDIjbLLV79247ubZt2xbgUQAAgHyT84TiZLz55pvmnnvuMQ899JDN0XnxxRfNq6++au66666Y7xk9erSpWbNmZFISMgAACK8iR1UgOWqWUn7N888/by688MLI/F69epktW7aYl19+ucR7fvnLX5qTTjrJjB07NjLvmWeeMf369TPbt2+3zVqJ1NwowNm6daupUaNGRo4NAAAES9dvVVIkcv3OWc1NxYoVTZs2bcycOXMi8/bv32+fd+jQwfc9O3fuLBHAlC9f3j7GitEqVapkC8E7AQCA8MpZzo2oG7hqak488UTTrl07O4bNjh07bO8p6dmzp2nYsKFtWpLzzjvP9rA6/vjj7Zg4K1asML///e/tfDfIAQAAZVtOg5vu3bubTZs2mREjRpj169eb4447zsyaNSuSZLx69epiNTW/+93vTFFRkX1cu3atqVOnjg1s/vCHP+TwKAAAQD7JWc5NIbTZAQCA/FAQOTcAAACZQHADAABCheAGAACECsENAAAIFYIbAAAQKgQ3AAAgVAhuAABAqBDcAACAUCG4AQAAoZJycLN3717zz3/+00yePNn88MMPdt4333xj784NAABQUPeWWrVqlTn77LPtvZ92795tzjzzTHPggQeaP/7xj/b5pEmTgt9TAACATNXc3HjjjfZO3t9//72pUqVKZH63bt3MnDlzUlklAABA7mpu3nnnHbNgwQJTsWLFYvObNWtm79YNAABQUDU3+/fvN/v27Ssx/+uvv7bNUwAAAAUV3Jx11llm/PjxkedFRUU2kXjkyJGma9euQe4fAABAUoocx3GSe8v/1dB06dLF6K3Lly+3+Td6rF27tnn77bdN3bp1Tb7atm2bqVmzptm6daupUaNGrncHAAAEfP1OKbhxu4JPnz7d/Pvf/7a1NieccILp0aNHsQTjfERwAwBA4cl4cKPamY4dO5oDDjigRMCjRONTTjnF5CuCGwAACk8y1++Ucm5OP/10891335WYrw3qNQAAgFxJKbhRZY+SiKN9++23plq1akHsFwAAQObHubnooovsowKbX//616ZSpUqR19Q1/MMPP7TNVQAAAAUR3Kity6250Xg23uRhDeh30kknmb59+wa/lwAAAJkIbp544onISMSDBw+mCQoAAOSdlLuCFyp6SwEAEO7rd0r3lpLnn3/ezJgxw94ZfM+ePcVeW7ZsWaqrBQAAyH5vqT//+c+md+/epl69eub999837dq1M7Vq1TJffvmlOeecc9LbIwAAgGwHNw899JB55JFHzIQJE2wi8ZAhQ8zs2bPNDTfcYKuLAAAACiq4UVOU2+VbPaZ++OEH+/+rr77aPPvss8HuIQAAQKaDm/r160dGKG7SpIlZtGiR/f/KlSttN3EAAICCCm46depkXnnlFft/5d7cfPPN5swzzzTdu3c33bp1C3ofAQAAMtsVfP/+/XZyb5w5bdo0e8PMli1bmv79+9s8nHxFV3AAAApPxu8KXsgIbgAAKDxZGedm165d9l5SGzdutLU4Xueff36qqwUAAEhLSsHNrFmzTM+ePc3mzZtLvKabauommgAAAAWTUHz99debSy+91Kxbty6Sf+NOBDYAAKDggpsNGzaYQYMG2RGKAQAACj64ueSSS8ybb74Z/N4AAACkKaXeUjt37rTNUnXq1DHHHnusqVChQrHXdRuGfEVvKQAACk/Ge0vpFgtvvPGGqVy5sq3BURKxS//P5+AGAACEW0rBzfDhw80dd9xhhg4dasqVS6llCwAAICNSikz27Nljb7VAYAMAAEIR3PTq1ctMnz49+L0BAADIRbOUxrK59957zeuvv25atWpVIqF43Lhx6e4XAABA9oKb//znP+b444+3///oo4+KveZNLgYAACiI4GbevHnB7wkAAEAA6OoEAADKZs3NRRddZKZOnWoHztH/43nxxReD2DcAAIDMBTcaFdDNp1GAQ24NAAAIze0XChm3XwAAINzX75Rybjp16mS2bNniu2G9BgAAkCspBTe6n5RGKY62a9cu88477wSxXwAAAJnvCv7hhx9G/v/JJ5+Y9evXFxvYb9asWaZhw4Z8FAAAoDCCm+OOO84mEmvya36qUqWKmTBhQtI7MXHiRDN27FgbLLVu3dquo127dr7Lnnbaaeatt94qMb9r167m1VdfTXrbAACgDAc3K1euNMo/PvTQQ82SJUtMnTp1Iq9VrFjR1K1b15QvXz6pHdA9qgYNGmQmTZpk2rdvb8aPH2+6dOliPv/8c7s+v27m3iaxb7/91gZEl156aVLbBQAA4ZTz3lIKaNq2bWsefPBB+3z//v2mcePG5vrrrzdDhw4t9f0KhkaMGGHWrVtnqlWrVury9JYCAKDwZLy31JNPPlmsCWjIkCHmoIMOMh07djSrVq1KeD2qgVm6dKnp3Lnz/3aoXDn7fOHChQmtY8qUKebyyy+PGdjs3r3bFoh3AgAA4ZVScHPPPffY/BpREKJaF90lvHbt2ubmm29OeD2bN2+2icj16tUrNl/PvcnKsahpTDfuvPbaa2MuM3r0aBvpuZNqhQAAQHilFNysWbPGtGjRwv5/5syZ5pJLLjH9+vWzgUQ2u4Kr1ubYY4+NmXwsw4YNs1VY7qR9BwAA4ZVScFO9enWbyCtvvPGGOfPMM+3/K1eubH788ceE16OaHiUgb9iwodh8Pa9fv37c9+7YscNMmzbN9OnTJ+5ylSpVsm1z3gkAAIRXSsGNghk1BWn64osvbDds+fjjj03Tpk0TXo96WLVp08bMmTMnMk8JxXreoUOHuO997rnnbD7NVVddlcohAACAkEopuNG4NEoeVs6MumbXqlXLzldy8JVXXpnUutQN/NFHH7VJyp9++qn5zW9+Y2tlevfubV/v2bOnbVrya5K68MILI9sGAABIepwbl3pGaVyZyZMnm1GjRpljjjnGjkx82GGH2TFwktG9e3ezadMm251bScQaKFAjHbtJxqtXr7Y9qLw0Bs78+fNtkxgAAEDa49y88MIL5uqrrzY9evQwTz/9tL0Vg4Ia9Zp67bXX7JSvGOcGAIDCk/Fxbu6++247orCakypUqBCZf/LJJ5tly5alskoAAIBApBTcqFnolFNOKTFfEdWWLVuC2C8AAIDsBTfqpr1ixYoS85UHk2zODQAAQM6Dm759+5obb7zRLF682N4h/JtvvjF/+ctfzODBg21vJwAAgILqLaUbWmo8mjPOOMPs3LnTNlFpsDwFN7rhJQAAQEHeFVw3vlTz1Pbt281RRx1lRy7Od/SWAgCg8CRz/U6p5sY7wrCCGgAAgILOuQEAAMhXBDcAACBUCG4AAECoENwAAIBQSSuhGP6+/vprs3z5ctt7TD3JWrZsaRo1akRxAQCQBQQ3AZsyZYrp16+fHQfIpbuaP/LII6ZPnz5Bbw4AAAQ5zk0hyuQ4N6qxadq0abHAxlW+fHnz1VdfUYMDAEA+3hUc/tQU5RfYyL59+3zvxwUAAIJFcBMg5daoCcqPam5atGgR5OYAAIAPgpsAKWlYuTUKZLz0fPLkyTRJAQCQBeTcZIByb9QEVa1aNbNjxw5bY0NvKQAACuDeUvCnQIZgBgCA3KBZCgAAhArBDQAACBWCGwAAECoENwAAIFQIbpBWr7B58+bZRwAA8gXBDVK+h5ZuNdGpUyf7qOcAAOQDxrlBIPfQ4t5ZAIBM4t5SyPo9tHTvrIULF1LyAICco1kKSXvvvfd8519++eU0TwEAco5mKaTdJOWlG4cuWrTItG3blpIFAASGZilktUnKS6+ddNJJ1OAAAHKGZikkpWXLlrZ2Jh4FOP3796eLOAAgJwhukBTdEPTqq68udTklGOvO6LnEODwAUDYR3CDpgOHpp58udTl1DW/RokXOSpdxeACg7CK4QaA5N64xY8bYWp5cBWD9+vWL7GeqzWTU/ABAYSK4QSDdwKOpR1WuPPDAA77j8CTTTEbNDwAULrqCI+1u4EVFRcZxnGLzlHT8yCOPmC5dutjaHiUiB1mTo33xW6/mN2nSxHd/Vq1aldA+MAIzAOQfuoIjq01St9xyS4keVFqub9++Gbn/VLxaFe1jdGDj7o9qdPyaprzNT5p0POnW/AAAcodmqSwIS+6GXzdwJQ7feOON5tlnny2xvIKMdPNeEsmn0fN333231GazP/3pTyWCIW+gpBqfxo0bmxkzZqTVJAcAyC2CmwwLU+7G66+/XqxWRIHO5MmTbVNPx44dSx3/JojaD7/aIz1v3769GTBggBk6dGjc93uDoehAya/Gx0vrLvQAFQDKAoKbAui1k0/HEh0AKKdGFOAoxyZegBNE9/BYtSfaLwVaifTkckdR7tOnT0LLu2iaAoDCQHCTg7tnp1N7kasmrlg1Jt5jUbAwbNiwuDUfiSb0+h2jnpdWM5Mo7fsbb7yR9PumTZtWkMEpAJQlBDc5yFFJtfYil01ciRyLLvqjR4+OuQ69Vto+xzrGWIm+2abaoUJvXgSA0HPKmK1bt6pdxT5mw2OPPeaUL1/eblOPep6KNWvWOOXKlbPrcSetT/OzQftdVFQU2bb2JfpY5s6dW2z//KZ4+xzrGMeOHVts28lM0esLatJ6lyxZkpWyBwA4SV2/qbnJMDXVfPXVV7aZRY96ni9NXEHl2yRzU814+xzrGIcMGVJqsm8sgwYNKnWfUsHdzwEgfxHcZIHyTE477bS0BrELuokr6HybZJKKq1Wr5ptTEy9ZOBVuN/VFixYlHOBoQEJN7vsvueSSUntekYMDAPmF4KZAuIGDLriiR7cbdqb5BR2xAivVTMUKJjSvW7dutqeSO67MrbfeGhk8L9lk4R49esQNWm6++WZbPm3bti016FIvNgVcq1evtpNb0/bb3/427j64gwMCAPIHt18oMAoCVGOiwCIbgU2sWy6MHTvWDB48OOb7lHCrgEHNSm5NSKwaGAUdaj7SIHuJcm+noLF3vN3to1/3lpHGtmnXrl3MdfmVZ6zjjw70FAjl6kahAFAWbNu2zdSsWdNs3brV1KhRI+6y1NyUwSauIG65cOKJJyaUa6TRfv3uPeWl9ScT2Mgf//hHWwbajgITBVreWi3V1ESX0fbt233XpcAqVnkm0tTG+DcAkF8IbgpILsa4SSfXR4FB7dq1M9J92xtcaTuqSSotcdvvWPRceTnxeAOoWEEOt2YAgPxBcFMgcjXGTbxbLiQikR5UyYoVXJVWq+WXt+RXwxPrvQqgFOSouS3abbfdZmupSC4GgNwj56YA+OV9ZCPPw2+78fJTYlEgdu211wayT24idapd6oPIW1LtkIJMPyofBUze/dP21LynQC/W9hJZBgDKsm2FlHMzceJE06xZM1O5cmV788MlS5bEXX7Lli3muuuuMw0aNDCVKlUyP//5z81rr71mwixXY9wk2gW8NBoPJ9naGzcJ2Us3xkxnrKCg8pbi1Ua59w9T8rKCIPdO5PFq3MJ0c1UAyAu5HPVw2rRpTsWKFZ3HH3/c+fjjj52+ffs6Bx10kLNhwwbf5Xfv3u2ceOKJTteuXZ358+c7K1eudN58803ngw8+yNsRioOQq9GJg9puIiMXe6dYoxHPmzfPyRf33ntvSiMjR5dfrkeeBoBCUTAjFI8bN8707dvX9O7d2xx11FFm0qRJpmrVqubxxx/3XV7zv/vuOzNz5kxz8skn2xqfU0891bRu3dqEPflXPXrcXBHVGmgMl0xTzYZ6Jbm1FKmOrZNs3s3tt9+eswELE1Vab7FYSdTRNW533313zkaeBoCwyllws2fPHrN06VLTuXPn/+1MuXL2+cKFC33f88orr5gOHTrYZql69eqZY445xtxzzz32YpBPgmxmcNel5g1dBN1bHrjNHZlswtC6lSir7aqZSDe+TKVJyC+RVwGMX9OT6BzI1YCFyQRsqdAxuyM0/+53v7PH5bdMPgVyAFBwnBxZu3atrV5asGBBsfm33nqr065dO9/3HH744U6lSpWca665xnnvvfdss9bBBx/sjBo1KuZ2du3aZauw3EnV/ZlslgqymcFvXcnciDJfjsO7TjUtuevwa9rxbiN6+XzinkeZuplnPh4zAORSwTRLJUs1CHXr1rW/6tu0aWO6d+9uhg8fbpuzYlFtg7Kr3alx48YFk/wbawC9INadyrbT3VZ0Iq9uvaDu1bGavbI9YGGy5ZOqRMb94ZYOAJC6nAU3GtxNF7MNGzYUm6/n9evX932Pekipd5TbXCFHHnmkWb9+vW3m8jNs2DDbbcyd1qxZYzIpyBtcJpKrkqlclGzdqFMD46lrebp3Tc+2TIzfE52Pxpg5AFBgwU3FihVt7cucOXOK/aLVc+XV+FESsWoOvL98v/jiCxv0aH1+1F1c/eG9U6Hc4NJdV3RuSroJvtkYvC8Z+VxDk+jnrPKJlUOUCm7ICQBpcHJIOTPKoZk6darzySefOP369bNdwdevX29fv/rqq52hQ4dGll+9erVz4IEHOgMHDnQ+//xz5+9//7tTt25d5+677867ruBB5Yv45b7o+YwZMzKWlxFrm+SBxP+cH3vsMZsz5JbX7bffbj+nWF3bE8nNocwBIPnr9wEmh5Qzs2nTJjNixAjbtHTccceZWbNm2Z5Qsnr16mJV/8qXUY2CukG3atXKNGzY0N4XSD168vGXfRA1EbEG0qtTp07GajriDd5XSLUr2f6c1aSm3mzRox9rYMrSbgyqruXR96dya2+UlwQASBy3X8hzubj1Qq5u9xBWfuXppQB+0aJFdoTu6LunU+4AUIC3X0D2cnjyMd+mLH6Gys1x83Pcm3e2bdvW3HLLLSXeqx5qzz33XMElF/sNYpmLu9oDKJuouclz7g0Vq1evbnbs2JHyzR6zfbNM+Jet22Ql0c1X8Wp4/G7Ima/UBOcd/FGjXOtxyJAhNmjW/x999NGCOBYAhVlzQ3CTxzRCcL9+/exFIlsXt1h3vNZ89WhCZim/RkGAn0IIMjXq8h/+8IdSl1OAo5y6fD4WAPmFZqkQ0K94N7Dx3m3ardLPVBV/tsa3QfL3rMr37uGJBjaiGpxYt1kBgHSRc5On4o0QHOS9qzJ1s0ykHlzGGy/n/vvvz8ucFdU4JRrYeO8VBwCZQHCTp/xqUPRcN12MV6OTLzfLRGoURPolFufbHcO9NYeaUhmO4Zlnnim1izwApILgJk/5jU6sqvwZM2YEfs+nWE1h2p5uX5GPNQVhprGb4t3aIXo8nGxz70jv1hz26NGjRBf2RCm/KNNNrQDKHoKbPKYB4aKDG91zKFM5MZm4WSbS7zoeTbUk2Q4A3MDjhhtusDc89dYcvv32277vSeR2FDqnlUfkDZiaNGnCwIUA0kJwk2Hp/BqNNVLwoEGDMjLuDcnE+UNNgRo0UcFsrhOLvYHHhAkTEnqPmkrV3dt7762uXbvGXL83YFLAoxodJSgLNToAkuaUMdm6t5ToXkPuPZr0qOfp3uNJ9y7S/KDuXRWtV69exban58gdfb5+96ZK9r5TWnbu3LkJvcdddsmSJU7//v2TvieW9tfdTvR5muz6OnXqlNbfEICyef0muMmQeIFJMrw3Y9RjJr/cg9pnBGvw4MG+F37NTyRIGTt2bIkAIVaw41021WnAgAFx9ynVG4lyM1GgbNtKcBNM4aRDFw6/L2f9ik1WpmppMrnPCE6sgCBe4OmtNfSrWXFf0//vvfde+x49phPUuFNp52msYC3RSbU/KJu8QXkytZHZ3C9kDsFNQIWT61qQbP/BUHOTv2IFBG7gGf2ln2zty/Dhw9OqUXEn1fyUJt3am0S3kw1l6aIW71g1b/r06XZSc6b3/9Hv8S6bTLlF1yq655A3QI+3n5n4rLQuBdux9gXBIrgJqHDSlU6TUrr5OqnyfoFkuhkMidNFwu8ir/nec0Vfrq1atQqkBibepO0pIPKeo8l8qadbS+TN60lEJi5s3r+VQr+oefOs/AKT6HPMe6wqh9KCVb1+++23++ZcaX5ptTE610o7J3r06FHsnHT30y8A8W4zWSqX++67z7n++uvj/lhItNy9gV5ZCpZTQXATUOEEIZUmpURrUIL+Q4j3BYbcitVkqC/tdHNkkp28QW86TabRuUCXXHJJ5MdAkM1TmQhCYl1sE72o5ZN4TZilBRPxLvCpBq3RgYm2k8nz+corr3QeeuihhGqSLr300oTXe8YZZzh/+9vfIkGj97vaDbhKO/5oZT342UrOTTCFkyuJ5L4EXbNDk1R+C6IpJ4hpxowZgVfre4MjPWobiV5s4zVPxbswptOsVVotQlBNZkFcyNyaAffi7V5kvY/5cF7l0xTdxOWW39lnn532ulXWXbt2TXj5bt26RT47b62UG5Ql2rS3JiRBEcFNQIWTK6UFGpkIREgmzn/pJuKmM2W7G3Z0k+4pp5wS82Lhl9ehAKO0Y4qVExJPIs0j0V3h4zW5JNJrzVuLkcy+JtJcxBQ7qAgqDy0bU6xAZ01Uk5ymfv36xT2HkjlPvc2Z2Uj0JrgJqHByKV6+TtCBiPvrhG7g4a290fuSea++2HU+6UsrGz31Yh2vu+1Ejl2vKwBMpalE7/U2tUXXbugxmfVqP0pLgI0OYPSeeGMLJZO0mkgQxhTOMtD589BDD5XanBcdEPk1lSkQis5jcvOV9Pfg9zeZyeRqgpuACidIqUSzsfIZgqy5iW5vd09MkonDU3tz//33R5JC3YBZn/npp59eUHkjQXVVjzfpCzvbOUypTrE+p7AHNvqM1HT58MMP503NylVXXRV4/lE2p/bt22ds3UH2bCS4CahwgpKJnk9BDO4Xq8swycThqb2JDnqjA+bopN58TyLPdHJpKtPRRx+dNwFOUIFNKkGDtu3mTGlya/2if/mrdkCvabDHZAPJ6PPT+z3oTlqn8lpiHUOQAVH030wizaFlbSqX5Gjq8RDcBFQ4Qchkom66g/vFat4Kch+RvdobXfj1WrJBb7rnUTapGj3XX9beSRfuXO+TG+AEEdioBsJtBowOTNwLuV7zBiaJBMXxaqHd5k83KIpVA6Jl4q07ugnVL9DyHpu7TdUAJZMs7DbNxPqb0TyVYyJBlHJ73P1z97WQcn1MglNQA8ES3ARUOEHI10TdWHk2+bSPSDxg9v46KqRgJVk6plS+XHXxDbpZyw0qUt2nICfdgyud98errS0tMMnEeeZNhs5WMnu8bt4KQhQEJdNbMDroSmZcKDcwiw6S3MDKfS3X551JYAryhzLBTUCFE4R87GKdyLgWud5H5M89xwo578itiQjy3lnuerOdD5TqpJoQ70XW++jWFuTj33ougnSVyXnnnZexptpUxz2L9Tm5tUSlnQP9+/fPSY1Q0N9NBDcBFU4YL0Sx8mw0kUxcWMJcQ5PqOexO8S5KbrlpgLVUv+z9yjxWs1AivbxiHYN3YMN4SeCxpnxNDs93hfa3Fd3zrkePHr41TW6gpNdiNZ2579Xr3iBPeUzRf3d6/8CBA0sEy5nqZZnM9btI/5gyZNu2baZmzZpm69atpkaNGlnb7tdff21WrFhhWrRoYRo1amRyZd68eaZTp04xXy8qKjJ//OMfza233prV/QKSMWXKFNO/f3+zb98+U758eTN69GjTtm1bU61aNbNjx46E/86863GVK1fODB061OzZs8fcf//9JV575JFHTJ8+fXzX96c//cncdtttZv/+/XbZMWPGmCuuuCLytz9t2rSYr7v7Hn0M0d8dV111lfnLX/5S6rENHz7c3H333QmXKQpbKteYr7/+2ixcuNB8++23platWqZDhw7F3hu9Tvd5sn9nubh+E9xkmU6O5cuXm5YtW8Y9KYJezrt806ZN7ZdrLLpYfPXVVzkNwoBs/WCI94Wdypd5afuV7n7r/Y0bN467DIENTFmvnHDKmFwO4pdol/Cgl4v3vlgTycRA/oqX40NTFMKKZqmgIr8A+dWY+NWQBL1crH1ZsGCBbaJSFXt0LQ41N0D+UxPYkCFD9APVPqdJGWG3LYnr9wFZ26syTk1H0UGE2vJVPe0NRoJezi/HoF+/fva9+jLU//X46KOPRvIXJk+eTJMUkOcGDx5sLr/8cpszIdH5EkBZRs5NGaq5iZVv4yY2KiEz1wnPAACkW3NTLu6rCIwCBjUBKQCRWDUkQS/n5VfbI5o3bNgwAhsAQChQc5NlifaUCHq5RHpKKQfntNNOS+JoAADIDrqCB1Q4mZZsN+4geHNuvEgiBgDkM5qlCoCCDNWiaEA9Pep5NnTp0sU8++yzZsCAAUk1aQEAUCholsqBdLpxB1VrQxIxAKCQUHOT5+J14041WFK+jB7jLeNtjiKJGAAQVvSWygHl2KjmxEs1N0oKzlTzVtABFQAA+YrgJgdS6cadaG2MbgLoV4MTZEAFAEA+I7jJEd1VWDk2ak7So/cuw4k0MyVbGxNUQAUAQL7j9gs5pMAiOriITvpVQOINfPxqY6ITk2PVxqin1F//+ld7uwWGagcAhBU1N3kkmWamZGtj3Nyc7t272/vRvP766xk+GgAAcoOu4HlETVFKDE525ODSRinOVddzAACCwl3BC1SyzUzxmre8Ur2DOAAAhYhmqTySqaRfekoBAMoSgpsC6kUVT7weVvSUAgCUJeTchECsHlbRN+ZM5g7iAADkE3JuypBYPay+//57c9ttt9nn6vp9yy23mMsuu8w4jpPrXQYAIKOouSlwsXpYRScmR78Wb/wcAADyDTfOLEP8koXjBTaJjJ8DAEAhI6G4wPklC48ZM6ZEwBONm2YCAMKK2y+EgJqXdGsFb7LwwQcfXCwXJxo3zQQAhBU1NyGhgEajGLu9oBTwrFq1ygwePDhSq+PippkAgDDLi+Bm4sSJplmzZqZy5cqmffv2ZsmSJTGXnTp1qu394530PpSkQGfs2LGRcXNUrsmOnwMAQKHJebPU9OnTzaBBg8ykSZNsYDN+/HjbxPL555+bunXr+r6nRo0a9nWXApx8Ez3GjHde9erVzfbt2+3jypUr7WvNmzePzNOj932pbM9vW6msFwCAQpPz4GbcuHGmb9++pnfv3va5gpxXX33VPP7442bo0KG+71EwU79+/bS2u3PPXnPAnr0mE5588ilz/fUDI4PqTZjwoJ3vzkuE+75evXomvb0rrrjCPPvss77bSma9AADkC123C2Kcmz179piqVaua559/3lx44YWR+b169TJbtmwxL7/8sm+z1LXXXmsaNmxoL94nnHCCueeee8zRRx/tu43du3fbydtPvnHjxqbxTTNMuUpVM3RkAAAgSPt37zRrxl9mtm7daltw8jbnZvPmzbZLcr169YrN1/P169f7vufwww+3tToKfJ555hkb4HTs2DHmmC2jR482NWvWjEwKbAAAQHjltObmm2++sTUwCxYsMB06dIjMHzJkiHnrrbfM4sWLS13HTz/9ZI488kjbFHPXXXclXHOzbtO3pUZ+qVi7dq3dH2+TkDvmTKJNUt5eTZ988okto2S2F8R6AQDIJ7p+N6hTK6Gam5zm3NSuXdteaDds2FBsvp4nmlNToUIFc/zxx9sxXvxUqlTJTtGqVjzATkFr2bypmTxxgh0BWLVSttv1w5Pta+68ROh9kx6eaNeX7PauuuoqW6vlt61E1wsAQD7Zm8Q1O+f3llIPqXbt2pkJEybY56qBaNKkiRk4cGDMhGIvXcCVb9O1a1ebnBzkvSnS4XcHbndetWrVzI4dO+yjumWLusK78/SY7J27o7fnt61U1gsAQD5I5vqd8+BGXcGVQDx58mQb5Kgr+IwZM8xnn31mc2969uxpm0+UOyN33nmnOemkk+xFWknHGsdl5syZZunSpeaoo47Km+AGAAAEJ5nrd867gnfv3t1s2rTJjBgxwiYRH3fccWbWrFmRJOPVq1cXu0/S999/b7uOa9mf/exnpk2bNjZnJ5HAJuyix7ZhTBsAQFmU85qbbAtrzc2UKVNK3EtKQaFuqsloxACAQldQzVLZFsbgRjU2TZs29e0xpQRi5fWQZwMAKCvX77y4txTSo6aoWF3BlXAdqycZAABhRHATAsqt8eYlRdfcKPkaAICyguAmBNTkpNwaBTJedoydyZNpkgIAlCnk3IRI9Ng2jGkDAAiLguoKjmBrcEgcBgCUdTRL5VnNy7x582LeBDRb6wAAoJAR3OTRODXqzt2pUyf7qOe5WAcAAIWOnJs88O6779pbSni7cyc7Po3fWDeMcQMACAvGuSkgql3RzUOjx6lJdnwav7FuGOMGAFAW0SyVQ6pt0S0T/AaJTnZ8Gr+xbhjjBgBQFhHc5OHIwgpSkh2fJnqsG8a4AQCUVeTc5JBfnowCm0WLFpm2bdumvE41ZzHGDQAgTMi5KRB+tS16nmpg467ztNNOY7wbAECZxSB+OaIaFjVLdenSxfaKorYFAIBgkHOTA9Hj0bz++uvUtgAAEBBybrKM8WgAAEgeOTd5jPFoAADILJqlsozxaAAAyCyCmyxjPBoAADKLnJscYTwaAAAyk3NDV/Ac1uAkMwIxAABIDM1SAAAgVAhuAABAqBDcAACAUCG4AQAAoUJwAwAAQoXgBgAAhArBDQAACBWCGwAAECoENwAAIFQIbgAAQKgQ3AAAgFApc/eWchwncgMuAABQGNzrtnsdj6fMBTc//PCDfWzcuHGudwUAAKRwHdfdweMpchIJgUJk//795ptvvjEHHnigKSoqCjyqVNC0Zs2aUm/HDsq5EHBOU85hwvlc2GWtcEWBzSGHHGLKlYufVVPmam5UII0aNcroNvRBEtxkHuWcPZQ15RwmnM+FW9al1di4SCgGAAChQnADAABCheAmQJUqVTIjR460j8gcyjl7KGvKOUw4n8tOWZe5hGIAABBu1NwAAIBQIbgBAAChQnADAABCheAGAACECsFNQCZOnGiaNWtmKleubNq3b2+WLFkS1KrLjLffftucd955dvRJjR49c+bMYq8r933EiBGmQYMGpkqVKqZz585m+fLlxZb57rvvTI8ePeygUQcddJDp06eP2b59e5aPJH+NHj3atG3b1o7QXbduXXPhhReazz//vNgyu3btMtddd52pVauWqV69urn44ovNhg0bii2zevVqc+6555qqVava9dx6661m7969WT6a/Pbwww+bVq1aRQYx69Chg/nHP/4ReZ1yzowxY8bY74+bbrqJsg7YqFGjbNl6pyOOOCI/z2n1lkJ6pk2b5lSsWNF5/PHHnY8//tjp27evc9BBBzkbNmygaJPw2muvOcOHD3defPFF9eBzXnrppWKvjxkzxqlZs6Yzc+ZM59///rdz/vnnO82bN3d+/PHHyDJnn32207p1a2fRokXOO++847Ro0cK54oor+Bz+vy5dujhPPPGE89FHHzkffPCB07VrV6dJkybO9u3bI2U0YMAAp3Hjxs6cOXOc9957zznppJOcjh07Rl7fu3evc8wxxzidO3d23n//ffu51a5d2xk2bBjl7PHKK684r776qvPFF184n3/+uXP77bc7FSpUsGVPOWfGkiVLnGbNmjmtWrVybrzxRs7pgI0cOdI5+uijnXXr1kWmTZs25eV3B8FNANq1a+dcd911kef79u1zDjnkEGf06NFBrL5Mig5u9u/f79SvX98ZO3ZsZN6WLVucSpUqOc8++6x9/sknn9j3vfvuu5Fl/vGPfzhFRUXO2rVrs3wEhWHjxo22zN56661ImeoC/Nxzz0WW+fTTT+0yCxcutM/1hVSuXDln/fr1kWUefvhhp0aNGs7u3btzcBSF42c/+5nz2GOPUc4Z8MMPPzgtW7Z0Zs+e7Zx66qmR4IZzOtjgRj8e/eRbOdMslaY9e/aYpUuX2iYS7/2r9HzhwoXprh7/38qVK8369euLlbPuMaImQLec9aimqBNPPDGyjJbX57F48WLK0sfWrVvt48EHH2wfdS7/9NNPxcpZ1c5NmjQpVs7HHnusqVevXmSZLl262Bvlffzxx5Szj3379plp06aZHTt22OYpyjl4ag5Rc4f33OWcDp5SAZQ6cOihh9oUADUzSb6d02XuxplB27x5s/3i8n5YouefffZZzvYrbBTYiF85u6/pUW24XgcccIC9cLvL4H/2799v8xJOPvlkc8wxx0TKsGLFijZIjFfOfp+D93PC//nPf/5jgxnlIigH4aWXXjJHHXWU+eCDDyjnAClwXLZsmXn33XdLvMY5HRz9mJw6dao5/PDDzbp168wdd9xhfvnLX5qPPvoo78qZ4AYow7909aU0f/78XO9KaOkioEBGNWTPP/+86dWrl3nrrbdyvVuhsmbNGnPjjTea2bNn2w4dyJxzzjkn8n8lyyvYadq0qZkxY4bt5JFPaJZKU+3atU358uVLZITref369dNdPf4/tyzjlbMeN27cWOx1ZeGrBxWfRXEDBw40f//73828efNMo0aNIvNVTmpq3bJlS9xy9vscvJ8T/o9+ybZo0cK0adPG9lRr3bq1eeCBByjnAKk5RH/3J5xwgq2p1aQA8s9//rP9v2oGOKczQ7U0P//5z82KFSvy7pwmuAngy0tfXHPmzClW3a/nqo5GMJo3b25Pfm85q51WuTRuOetRf1j6snPNnTvXfh76hYH/606vwEbNIyoblauXzuUKFSoUK2d1FVe7urec1dziDST1q1ndndXkgth0Lu7evZtyDtAZZ5xhz0fVkLmT8u6UD+L+n3M6MzTMxn//+187PEfefXcEmp5chruCq9fO1KlTbY+dfv362a7g3oxwJNbbQd0DNenUHDdunP3/qlWrIl3BVa4vv/yy8+GHHzoXXHCBb1fw448/3lm8eLEzf/5823uCruD/85vf/MZ2p3/zzTeLdefcuXNnse6c6h4+d+5c252zQ4cOdoruznnWWWfZ7uSzZs1y6tSpQ1fwKEOHDrW90FauXGnPVz1Xz7033niDcs4wb28pzung3HLLLfa7Q+f0v/71L9ulW1251esy3747CG4CMmHCBPuharwbdQ3XOCtIzrx582xQEz316tUr0h3897//vVOvXj0bTJ5xxhl2/BCvb7/91gYz1atXt90Le/fubYMm/B+/8tWksW9cChZ/+9vf2m7LVatWdbp162YDIK+vvvrKOeecc5wqVarYLzd96f30008Us8c111zjNG3a1H4n6Atc56sb2FDO2Q1uOKeD0b17d6dBgwb2nG7YsKF9vmLFirws5yL9E2xdEAAAQO6QcwMAAEKF4AYAAIQKwQ0AAAgVghsAABAqBDcAACBUCG4AAECoENwAAIBQIbgBAAChQnADAABCheAGAACEygG53gEASNdpp51mWrVqZSpXrmwee+wxU7FiRTNgwAAzatQoChcog6i5ARAKTz75pKlWrZpZvHixuffee82dd95pZs+enevdApAD3DgTQChqbvbt22feeeedyLx27dqZTp06mTFjxuR03wBkHzU3AEJBzVJeDRo0MBs3bszZ/gDIHYIbAKFQoUKFYs+LiorM/v37c7Y/AHKH4AYAAIQKwQ0AAAgVghsAABAq9JYCAAChQs0NAAAIFYIbAAAQKgQ3AAAgVAhuAABAqBDcAACAUCG4AQAAoUJwAwAAQoXgBgAAhArBDQAACBWCGwAAECoENwAAIFQIbgAAgAmT/wdcHdNkrUHvswAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final estimate: 0.558\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.351385\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
}