{
 "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": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "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
}