ian-whitestone/simulations.ipynb

## simulations.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T21:09:06.058959Z",
     "iopub.status.busy": "2021-01-31T21:09:06.058623Z",
     "iopub.status.idle": "2021-01-31T21:09:08.381710Z",
     "shell.execute_reply": "2021-01-31T21:09:08.381070Z",
     "shell.execute_reply.started": "2021-01-31T21:09:06.058904Z"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as ticker\n",
    "\n",
    "import numpy as np\n",
    "from scipy.stats import beta\n",
    "import pandas as pd\n",
    "import pymc3 as pm\n",
    "\n",
    "from dask.distributed import Client, LocalCluster\n",
    "import dask.delayed\n",
    "\n",
    "import logging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Session level randomization with independent sessions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T21:09:08.383547Z",
     "iopub.status.busy": "2021-01-31T21:09:08.383327Z",
     "iopub.status.idle": "2021-01-31T21:09:08.390139Z",
     "shell.execute_reply": "2021-01-31T21:09:08.389538Z",
     "shell.execute_reply.started": "2021-01-31T21:09:08.383522Z"
    }
   },
   "outputs": [],
   "source": [
    "num_users = 10000\n",
    "\n",
    "sessions_per_user = np.random.geometric(0.5, size=num_users)\n",
    "\n",
    "# each user has some baseline conversion rate, we'll say ~10%\n",
    "baseline_conversion_rates = np.random.beta(100, 900, size=num_users);\n",
    "\n",
    "# treatment group will get a +2% lift in conversion (20% relative increase)\n",
    "conversion_uplifts = np.random.beta(20, 980, size=num_users);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T21:09:08.393825Z",
     "iopub.status.busy": "2021-01-31T21:09:08.393481Z",
     "iopub.status.idle": "2021-01-31T21:09:08.932413Z",
     "shell.execute_reply": "2021-01-31T21:09:08.931929Z",
     "shell.execute_reply.started": "2021-01-31T21:09:08.393798Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8cAAAEICAYAAAByLoYSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAArIklEQVR4nO3dd7gkRb3/8c8HdllyxrACrgIiwZxQUPeKCZZkFjGs6crPi8oVL3cJKqDAGsCEXgMqkhaQJMkreBWMgIBLElDCwpIkLlFwhe/vj6rZ7TM78ZzpCaffr+c5z5mZ7umqTtX17aqucUQIAAAAAIAqW27QGQAAAAAAYNAIjgEAAAAAlUdwDAAAAACoPIJjAAAAAEDlERwDAAAAACqP4BgAAAAAUHkEx31g+9W2rxt0PiYj27Nt/67w/mHbz+7Rsve1fWR+PcN22J7So2VvmPO6fC+WByAps0yYLIplG4DJjfpGe7YX2H59fj2mfLT9FtsL8zZ8UQ/S+n2ny8n1zo0nmmbdMs+3/ZFeLrOw7MNs/78ylt1PlQ2ObW9j+w+2H7B9Xz5YX1ZGWhHx24jYtIxlY6yIWDUibmw1j+2Ztm/tYFmHRERPCpBiwZuXfUvO6xO9WD4wjPJx/49cqbjf9tm2N+hnHjopE6qml2XbRNWXjcAosf0e25fkMu4O2z+3vc2g81U0CvWNXjdATESD8vGrkvaIiFUl3T+RfNreUdJDEfHnXuS1bOMIpL8qaV/bK5SVp36oZHBse3VJZ0n6lqS1JT1D0oGSHh9kvoaRk4EcJ4MsJIehgAYmiR1zpeLpkv6uVO6iRMNSfg1LPoAy2P60pK9LOkTSUyVtKOk7knYeYLbGGJZzcMRbrZ8p6eoeLWt3Scf0aFlDJyLukHStpJ0GnZeJqGRwLOk5khQR8yLiiYj4R0ScGxFX1Gaw/SHb1+TWjl/Yfmb+3La/Zvsu2w/avtL2lnna9rb/Yvsh27fZ/kz+fExLpe3N8t2YRbavtr1TYdpRtr+dW1gesn2R7Y3apV0vL/9Q2xfneX9me+3C9K1yy/ki25fbnln33YNt/17So5KW6ZJY39Uj5/uL+fW6ts/Ky77P9m9rAbbt6bZPsX237Ztsf7KwjANsn2z7WNsPSprdIN11bJ+R1+liSRs1y1ej/WF7FUk/lzQ93+l9OOdpmbTzZ8fWZeFDtm/Pd4g/02j98/sl+9z2MUoXzTNzenu77i5pzsMZeXtdb/ujddvlJNtH53W52vZL67cNMMwi4jFJJ0vavPaZ7Vm2/5zP54W2DyhMWzGfj/fmsuRPtp+ap61h+4f5PLzN9hfdpPJVVyY0LV/z9OfaPi+fh9fZfmez9bG9tu0f5/LgftunF6Z9NJ/H9+Xzenpdfna3/be8Xt92Mi2/37Iw73pOLe9Pye93sD0/z/cH288vzLvA9n/bvkLSI7an5Pe35XW9zva2ed4xZZvtnXK5ssip/N+sbrmfsX2FU0+rE22v2GSbzHbqhfU12/dKOsD2RrZ/lffjPbaPs71mnn+ZsjF/3ur6NNv2jXmdbrK9W7N9BJTF9hqSDpL0HxFxakQ8EhGLI+LMiPivPM8021/PZcTt+fW0PG2m7Vtt7+VUp7vD9gfztFfYvrNYpjl17b0iv17O9hzbN+Tz6iTn+p2X1i0+bPsWSb/ysvWNhudQXu7+tm/OeTo6r2dxuR+wfUs+l/drsX2Osv0/ts+x/Yikf3OL8l7Sb/L/RbkseGVeTsO6eIP0lukR6LFdpWt1vBPzel9m+wVNlnWA07Vnmu2HJS0v6XLbNzTKp+2NbV+Qy8d7bJ/YZLkrSHqdpAsKn73c9h9zWXeH7SO8bKvr9nl/3WP7K15an26aru1XOV0zH8j/X9VqXQvvlxwrtg+W9GpJR+R1PSLP0+46eb6kWY3SGxkRUbk/SatLulfSTyRtJ2mtuuk7S7pe0maSpkjaX9If8rQ3SbpU0pqSnOd5ep52h6RX59drSXpxfj1T0q359dS87H0l1U6UhyRtmqcflfP28pz2cZJOaJd2g3U8X9JtkraUtIqkUyQdm6c9I6exvdINkjfk9+sVvnuLpC1yHqY2WH5I2rjw/ihJX8yvD5X03byuU5VOLue0LpX0ubzuz5Z0o6Q35e8dIGmxpF3yvCs1SPcESSflddoyr+PvGuWrk/1R+N4yaefPattsRl72vJz28yTdLen19evfKA1JC2rz1i1vSn7/G6U7zitKemFe9usKeXss76/l8/a9cNDnEX/8tfsrHveSVlYqc48uTJ+Zz6XlJD1fqWV5lzztY5LOzN9bXtJLJK2ep50m6Xv5XHyKpIslfSxPm92iTDhKzcvXVSQtlPTBPO1Fku6RtHmTdTtb0om5bJkq6bX589fl771Y0jSllvLf1OXnLKVyfMN8rr85T/uRpIML8/6HpP/Nr18k6S5Jr8jb4wN5+04rbOv5kjZQKr82zeszPU+fIWmj/PoALS3bniPpEaXrwFRJeytdo1YoLPdiSdOVelpdI2n3JttktqR/SfpE3oYrSdo4L3uapPWUyrqvNzpG8vum16e8jx7U0uvl0yVtMejjnL/q/Ul6cz7Wp7SY5yBJFyqVUetJ+oOkL+RpM/P3D8rn3fZKjRFr5ek3SHpDYVk/lTQnv/5UXu76+bz6nqR5edqMXMYcnc+XlQqfTWl1Dkn6UD73ny1pVUmnSjqmbrk/yMt8gVJvy82arPtRkh6QtHU+j1dU6/J+SR4Ly2haF2+Q3kwtW69bUrZoaR3v7Xl7f0bSTcr12wbzHltYTvEa0iif8yTtV1jPbZrkcQtJj9R99hJJW+X1m6FUvu5Zl/avlcreDSX9VdJHWqWb571f0vvycnfN79fJ088vLKN+XcesX3He/L7tdVLSWyVdNuhzdCJ/lWw5jogHJW2jpSf63U5395+aZ9ld0qERcU1E/Eupy8wL8x2rxZJWk/RcSc7z3JG/t1jS5rZXj4j7I+KyBslvpVTozI2If0bEr5QqSrsW5jktIi7OaR+nFCzVlt8s7UaOiYirIuIRSZ+V9M58J/K9ks6JiHMi4smIOE/SJUqFc81REXF1RPwrIha3SKORxUoF7jMj3Un9baQz5mVKAfhBed1vVNr+7y58948RcXrO1z+KC815f5ukz0W6S3uVUmW7VT7a7Y+ipmkXHJjTvlLSjzV2v42L0zOYW0v674h4LCLmSzpS0vsLs/0u768nlLrkNLzjCQyh020vUqoovUHSV2oTIuL8iLgyn3NXKF3sX5snL5a0jlKl5ImIuDQiHszl9PZKFYhHIuIuSV/T2HKklWbl6w6SFkTEj3O592elm4rvqF+A7acr3VjdPZctiyOi1hqwm6QfRcRlEfG4pH0kvdL2jMIi5kbEooi4RaniU8vD8XXr8Z78mST9u6TvRcRFeXv8RKlyulVh/m9GxMJcfj2hVHHe3PbUiFgQETc02B7vknR2RJyXy/qvKlV+iy0N34yI2yPiPqUbFi9cdjFL3B4R38rb8B8RcX1e9uMRcbekw7V0HzfS7vr0pKQtba8UEXdERK+6OwLdWEfSPbkcaWY3SQdFxF352D9QKWCpWZynL46IcyQ9rHRTS0pl4a6SZHs1peN/Xp62u6T9IuLWXMYcIOntHtuF+oBcPjaqyzQ7h3aTdHhE3BgRDyuVXe+uW+6B+by+XNLlal0X+VlE/D6fx4+1Ke8baVUXH49LI+LkXM4drhRQbtXmO51YrNT1enpez981mW9NpcawJfJ17cJcXi5QutFRv02+FBH35evF17W03tks3VmS/hYRx+TlzlPq6rzjhNYy6eQ6+VBe15FVyeBYkvLJNjsi1ldqgZyudNBJ6WD7Ru7msEjSfUotn8+IFMweIenbku6y/X2nZ5ilFLhtL+nm3NXhlQ2Sni5pYUQ8WfjsZqW75TV3Fl4/qhRMq03ajSysS2OqpHXz+r2jtn55HbdRCmgbfbdbX1G623du7goyJ3/+TKXuzMV091V6VqeTdNdTulNVv17NdLI/ijpZ5/q0pzebsQvTJd0XEcVCs90xsaKH5FkioI1dImJNpYrIHpIusP00aUn3wV87PWbxgFJlaN38vWMk/ULSCU5dEr9se6pSOTJV0h2FcuR7Sq0znWhYvublvqKufNpN0tMaLGMDpXP2/gbTpqtQLuVK5r3qoIxXCpRXzttlhlIQelohf3vV5W8DjS2DlpRPEXG9pD2VKs532T7Bhe7dLfL7ZF5OJ/ltZEw5avupOe3bnB5ZOVZL93EjTa9PkW70vkvpOLnDqXv8c1ssCyjLvZLWbXMdHnNuadk6w711wXXx3Dpe0ludumHXWuJqy3qmpNMK58c1SjfD2tal2pxDjfI7pW65EykLWpX3jTSti7f4TivF8vFJSbeqN3W4vXO+LnZ6POVDTea7X6mBawnbz3F6DPHOXD4eomW3SbN6Z7N06/dj7Xvj3W5FnVwnV5O0qAdpDUxlg+OiiLhWqQtI7VmvhUpd9NYs/K0UEX/I838zIl6i9OzccyT9V/78TxGxs1Il7XSl7r/1bpe0gccOcrWhUvfgTvLaMO0miqPCbqh0l+mevH7H1K3fKhExt5hUm6w8qtTdsWbJiRERD0XEXhHxbKWH8j/t9KzbQkk31aW7WkQUW6xbpXu3Ujek+vVqqMX+aJZGu3VWg7Rvz68fUZPt0cGyb5e0dr47XFx2R8cEMApya+epSpW42miux0s6Q9IGEbGG0uMYzvMvjogDI2JzpVbMHZR6UyxUajFdt1COrB4RW0wwiwslXVBXPq0aEY1+lmKh0jm7ZoNptytVICRJTuMcrKMOzudIPUNOUmoZ2FXSWYWbZguVulwX87dybhVYsoi65R0fEdvk/ISkL3WQXyuVc+Mtf+rLukPyZ8+LiNWVWobdYv6W16eI+EVEvEHpZu61Sr2PgH77o1I5tEuLecacWxpbZ2gpIv6iFNBsp7E9SKR0jmxXd46sGBHFc7ZpnaPFOdQov/9S6v48HvV5aFreN8lvy7p4nTF1sNzTcL26eTYoTF9OqVt6R/ujYJl8RsSdEfHRiJiu9DjQd9z455euT0m7GKT+j9I+2CSXj/tqbPk4Jt8qHEMt0q3fj7XvNSrTu627dnKd3EypV8HIqmRwnB8m38v2+vn9BkoVkQvzLN+VtI/tLfL0NWy/I79+Wb77NVXpoHpM0pO2V7C9m+01cpeNB5W6rtS7SCmw3Nv2VKeBRnZUepa2Xb4bpt3iK++1vbntlZWeazk5V76OlbSj7TfZXt5p4JuZte3RofmS3pO//2YVuoE4DRqzca5kPaBUGX5S6bm1h5wGiVkpf3dLd/gTWjnvpyoN8rKy7c2VnrtbRpv98XdJ6zgPNNGlz+a0t1B65qI2AMJ8pUET1s6tYnvWfe/vajCwWV6vhUrPIh2a98XzJX1YaT8Bk4KTnZWe0b0mf7yaUgvsY7ZfrlQJrM3/b7aflys5Dyrd3Hsy0qMk50o6zPbqToPIbGS7Vfe8Tpwl6Tm235fL5qm5zN2sfsach58rVUbWyvO+Jk+eJ+mDtl+YW30OkXRR7jLXieOVWnZ209gK8Q8k7Z6vAba9itMAN6s1WojtTW2/LufhMUn/UOPrxUmSZtneNl9b9lKq9DeqgI7HakrdRR/IlcL6G7r1ZWPT65NTK/TO+YbD43m5ra6BQCki4gGl8VO+bXuXXC+Yans721/Os82TtL/TwHrr5vm7ua4fr/R88WuUnjmu+a6kg710oNj1ctnaVptzaJ6k/7T9LNurKpVdJ0brruPdaFreKzV+PKmxZUHTungDf1XqUTcrl2P7Kz1WUvQS2291au3fU2n9L1R3lsmn7XcU6s/3KwWUy5RLEfFPSb/U2G7Tqyld3x52asFvdDP2v/J1ZgOl4+HENumeo3Qte4/TwFrvUmpQO6vBsudLeo3Tb2GvodSVvqi+fO7kOvlapevjyKpkcKzUH/4Vki5yGkXvQklXKVUKFBGnKd1hP8Gpm8NVSnfvpDSY1w+UDsSblbrW1J6he5+kBfk7uytVbsbIJ8eOeXn3KA3C9P7cet1Oq7QbOUapRfxOpS6Nn8x5WKg00MG+Sif6QqUKSzfHw6fyeixSWs/TC9M2USoAHla6u/qdiPh1Dm53UOoqeJPS+h8pqZsgdQ+lbjx35nX7cYt5G+6PvK3nSbrRqVtIN91qLlC6+/d/kr4aEefmz49RulO2QKniXj9a4aFKF8lFLoxyXbCr0kAItyt1o/x8RPyyi3wBw+pMpxE/H5R0sKQPxNJn3D4u6SDbDylVHIu9bZ6mNLr1g0rB9AVa+hMY71ca1O8vSuXhyRr7WEjXcgvtG5We+b1dqYz5kpatYNW8Tylgv1ZpoKw983J+qTTGwylKgwJupM6fh1ZEXKR083O6ChWMiLhE0keVHq25X6kcmt1iUdMkzVUqZ+9U6kFTX/FRRFyn1Jr7rTzvjko/v/XPTvPcxoFKg5M9oDSI2al108eUjW2uT8tJ+rTS/rlPqRLWqDIJlC4iDlM6HvfX0mN1Dy2tD31R6Xn5KyRdKemy/Fmnas/k/ioi7il8/g2lFthzc9l5oVKdthOtzqEfKZWxv1Gqoz2mNLherzQt7yPiUaXrw+9zWbBVm7r4GPlmxceV6pS3KZWht9bN9jOlG4+1wareGl2OqdMon0rj6VyUr3NnSPpUpDF1Gvmexj53/hmlmwQPKdXvG410/TOlwWznK5WhP8yfN0w3Iu5VqmvvpRQn7C1ph7pjqLY+5+U0r8hp1AfQ31B6nv1+299sd510Go9jc42NCUaOIzrpSYpRY/t8pRHojhx0XgAAAIBBcPrZqI0j4r1DkJffS9oj0mBWk4rtwyTdEBHfGXReJoIBfQAAAACgZBGx9aDzUJaI2GvQeeiFqnarBgAAAABgCbpVAwAAAAAqj5ZjAAAAAEDldfXM8brrrhszZswoKSsAqurSSy+9JyLqf5Nw0qDsBFCGyVx2Um4CKEO7crOr4HjGjBm65JJLJp4rACiwffOg81Amyk4AZZjMZSflJoAytCs36VYNAAAAAKg8gmMAAAAAQOURHAMAAAAAKo/gGAAAAABQeQTHAAAAAIDKIzgGAAAAAFQewTEAAAAAoPIIjgEAAAAAlUdwDAAAAACovCmDzkAvzJhzdtffWTB3Vgk5AQCMquK1hGsEACTtysbadMpNTAa0HAMAAAAVNmPO2eNqbAImG4JjAAAAAEDlTYpu1QAAAAA6R0sxsCxajgEAAAC0RfdrTHYExwAAAACAyiM4BgAAAABUHsExAAAAAKDyCI4BAAAAAJXHaNUAAIxTbWCaBXNnDTgnANBbDLyFKqLlGAAAAABQeQTHAAAAAIDKo1s1AAAAgAkpdsPmUROMKoJjAAC6wHN4ANAa4zFgVNGtGgAAAABQeQTHAAAAAIDKo1s1AAAdoDs1AACTGy3HAAAAAIDKIzgGAAAA0Dcz5pxNbxwMJbpVAwAqi8oZAACooeUYAAAAAFB5BMcAAAAAgMojOAYAAAAAVB7PHAMAAABgHAZUHi3HAAAAAIDKo+UYAIA+KbbKLJg7a4A5AQAA9QiOAQAAAHSsV92vuWGIYUO3agAAAABA5REcAwAAAAAqj27VAIDKYURWAABQj+AYAAAAQM/xTDFGDd2qAQAAAACVR8sxAAAAgFLxOAtGAcExAABNdFqZa9R1kO6EAACMFoJjAABK1ijIrn1G4AwAwHDgmWMAAAAAQOXRcgwAAABUBM/+As0RHAMAUIfKIwAA1UNwDAAAAExio3TDj8EMMUg8cwwAAAAAqDyCYwAAAABA5REcAwAAAAAqj+AYAAAAAFB5BMcAAAAAgMojOAYAAAAAVB4/5QQAAABgoEbp56YwedFyDAAAAACoPFqOAQAYEsWWkwVzZw0wJwAAVA/BMQAAA0RXQgAAhgPBMQAAPUSwCwC9Ra8a9AvPHAMAAAAAKo/gGAAAAABQeQTHAAAAAIDK45ljAAAAAEOHMRzQb7QcAwAAAAAqj5ZjAAAAYJKh1RXoHi3HAAAAAIDKo+UYAFAJtKIAAIBWaDkGAAAAAFQewTEAAAAAoPIIjgEAAAAAlUdwDAAAAACoPIJjAAAAAEDlERwDAAAAACqPn3ICAIy04k80LZg7a4A5AQAAo4yWYwAAAABA5REcAwAAAAAqj27VAAAAwAjj8RKgN2g5BgAAAABUHi3HAACMGFqJADRTLB8AdIfgGAAw6UyGymFtHQh+AQDoD4JjAMCkMRmCYgBAe/SgQRl45hgAAAAAUHkExwAAAACAyiM4BgAAAABUHsExAAAAAKDyCI4BAAAAAJVHcAwAAAAAqDyCYwAAAABA5fE7xwAAAABGAr9njzLRcgwAAAAAqDyCYwAAAABA5REcAwAAAAAqj2eOAQCYBIrP4S2YO2uAOQEAYDQRHAMAMMQYfAYAgP6gWzUAAAAAoPIIjgEAAAAAlUdwDAAAAACoPIJjAAAAAEDlMSAXAAAAMIIYsA/oLYJjAABGGJVjAAB6g+AYAAAAGBHcEGuO33vHRPHMMQAAAACg8mg5BgCMJFpPAABALxEcAwAAAEOOG4JA+ehWDQAAAACoPIJjAAAAAEDlERwDAAAAmFRmzDmbrujoGsExAAAAAKDyCI4BAAAAAJXHaNUAgKFX6xq3YO6sAecEAPqHbsFAf9FyDAAAAACoPIJjAAAAAEDlERwDADDJMEorAADd45ljAMDIIOADAABlITjuwngqZQweAwAAAADDj27VAAAAAIDKIzgGAAAAAFQewTEAAAAAoPIIjgEAAAAAlUdwDAAAAACoPEarBgBgkir+ygK/ngAAQGu0HAMAAAAAKo/gGAAAAABQeQTHAAAAAIDK45ljAAAAACOrOL4CMBEExwCAoURlBwAA9BPdqgEAqIAZc87mhgMAAC0QHAMAAAAAKo/gGAAAAABQeQTHAAAAAIDKY0AuAAAAYIgwPkDvFLflgrmzBpgTjAJajgEAAAAAlUfLMQAAAIBJr1GLPK3JKCI4BgAMFboTAgCAQaBbNQAAAACg8giOAQAAAACVR3AMAAAAAKg8njkGAAAABozxFoDBo+UYAAAAAFB5BMcAAAAAgMojOAYAAAAAVB7BMQAAAACg8giOAQAAAACVx2jVAICBqY3OumDurAHnpDqKI+Ky3QEAWIqWYwAAAABA5REcAwAAAAAqj+AYAAAAAFB5BMcAAAAAgMpjQC4AwMAVB4kCAAAYBIJjAAAAYEC4OTg8+AUF0K0aAAAAAFB5BMcAAAAAgMojOAYAAAAAVB7BMQAAAACg8giOAQAAAACVx2jVAABUFCOzAqg6RgtHES3HAAAAAIDKIzgGAAAAAFQewTEAAAAAoPIIjgEAAAAAlceAXAAAVFxxQBoG5wJQdZSJ1UVwDAAAAPQRIyQDw4lu1QAAAACAyqPleAiN524iXT4AAAAAYPxoOQYAAA3NmHM23T8BAJVBcAwAAAAAqDyCYwAAAABA5fHMMQCgr+imCwAYFbVrFuP7VAPBMQAAWIKbFwCAqqJbNQAAAACg8giOAQAAAKBDjOQ/edGtGgBQOioRAABg2BEcAwBKQ1AMAABGBd2qAQAAAACVR8sxAAAAUJJiDxp+DggYbrQcAwAAAAAqj5ZjAAAAAGiBMTSqgZZjAAAAAEDl0XIMAABaatRiwrOTAIDJhpZjAAAAAEDlERwDAAAAACqPbtUAAAAA0KV2P9NVm85jKKOD4BgA0FOM6Fk9VACBsSgHgdFEt2oAAAAAQOXRcgwAAAD0AS3Kkxc9aCYHgmMAANA1KvkAgMmGbtUAAAAAgMqj5RgAMG7tRuoEAAAYFQTHAAAAANADPHIy2giOAQA9QYUAAIBl0ctqdBAcV9h4KrKc0AAAAI1xkxAYbQzIBQAAAACoPFqOAQBdo3UEAABMNgTHAACgJxo9V8ezdgCAUUG3agAAAABA5dFyDABoqdbyR6sfutGq6z2tyQCqimvqcCM4BgAAAIAB4Ybh8CA4BgAAAMaJAQqByYPgGAAAAOgSQTEmguNnODEgFwAAAACg8mg5BgB0hLvcmCiOIQDAMCM4BgAAAIAhwiBdg0FwjNKNp6WAQgAAAABAPxEcAwCWQfdXAAD6j+vvYDEgFwAAAACg8mg5BgAswR1rDEqj5+t45g4A0E8ExwBQcQTEGKRGxx/HJIYVxyYwuREcA0CF0BIHAN0jKMYg1Y4/rtvlIzgGAAAjiZs9AKqKgLkcBMcAAGDoEQgDqLp2PRgImCeO4BiTBr+nDADVQBdX9APHGUZVu5uJBNHNERwDwCTXrIJHxQ+TXaNjnMog2qFsxCihNbm3CI6BLtFCDQAAgMmuij+xR3AMAABGHq196BTHCtC9qpw3BMfAkKKFGhNVlQsZUKZm5xHl7XCjSz0wVhm/KT8Zu2wTHAMVRxA+fNpV6hpdjAiEgc60Or/Gcx5NZOCbyd49sd8oB4HBaHfuteqS3U052I+bXo6Izme275Z0c4/SXlfSPT1aFumQDumMdjqbRsRqfUhnICZYdvZrH3RimPIikZ92hik/w5QXafLk55kRsV6vMzMMelznnKhhO15qhjVf0vDmjXx1ZzLmq2W52VVw3Eu2L4mIl5IO6ZAO6fQrnVE0TNtmmPIikZ92hik/w5QXifygO8O6f4Y1X9Lw5o18daeK+VqujIUCAAAAADBKCI4BAAAAAJU3yOD4+6RDOqRDOn1OZxQN07YZprxI5KedYcrPMOVFIj/ozrDun2HNlzS8eSNf3alcvgb2zDEAAAAAAMOCbtUAAAAAgMojOAYAAAAAVF7fg2Pbb7Z9ne3rbc8pMZ0f2b7L9lUlprGB7V/b/ovtq21/qqR0VrR9se3LczoHlpFOIb3lbf/Z9lklprHA9pW259u+pMR01rR9su1rbV9j+5UlpLFpXo/a34O29+x1Ojmt/8zHwFW259lesaR0PpXTuLrX69Lo3LS9tu3zbP8t/1+rl2kOk06OFyffzOXkFbZfXPjupfmzV+bPptj+pe2VJ5CnlseV7Wm2T8z5ucj2jPz51jkvl9jeJH+2pu1zbY/r+tLu2Ct720zk+LT9gTzP32x/IH82zfb/5nX6eGHe79fyPo78vCNvnydtN/0pCze53to+Lm+nQwqf7W97l3Hm5yu5jL3C9mm21+xHfprk5Qt5WfPzcTi9yXf7sq8K0/ayHbbX7Vd+MFar/ZOnNyxb8rRS908HeZtp+wEvvW58rjCt5+d5Yf629dxBbLcO89X3beYO6utucj3N0/bJn19n+035s/Vs/y5vr10K8/7MTcq3ceZrtu27C9vrI4VpZR//TeOOQWwvRUTf/iQtL+kGSc+WtIKkyyVtXlJar5H0YklXlbg+T5f04vx6NUl/LWN9JFnSqvn1VEkXSdqqxPX6tKTjJZ1VYhoLJK1b1vIL6fxE0kfy6xUkrVlyestLulPpB8Z7vexnSLpJ0kr5/UmSZpeQzpaSrpK0sqQpkn4paeMeLn+Zc1PSlyXNya/nSPpS2cfGMPw1O14kbS/p5/nc30rSRfnzwyVtI2l9Safkzz4xkeOgk+NK0sclfTe/frekE/PrU3NetpF0WP7sq5JmlnXslb1txnt8Slpb0o35/1r59VqSdpK0v9LN6D/meV8g6YcTyM9mkjaVdL6kl7Y4tpa53kp6vqQj8zznSVpD6Vp25gTy80ZJU/LrLzXZPj3PT5O8rF54/cnacTuofZU/30DSLyTdrAbXvbLyw19n+6cwvVnZUvr+6SBvM9WgTlbWeV5Yftt67iC2W4f56vs2Uwf1dTW/nm6e8zJN0rNyHpdXKsfeq3RdPD/Pu6OkA3qcr9mSjmjw3X4c/03jjkFsr363HL9c0vURcWNE/FPSCZJ2LiOhiPiNpPvKWHYhjTsi4rL8+iFJ1yhVNHudTkTEw/nt1PxXykhqtteXNEvSkWUsv59sr6F0wfmhJEXEPyNiUcnJbivphoi4uaTlT5G0ku0pSif+7SWksZnSxe3RiPiXpAskvbVXC29ybu6sdCND+f8uvUpvyDU7XnaWdHQ+9y+UtKbtp0tarLTfV5a02KmFbkdJR08wH+2Oq+L+OVnStrbdID8bSdogIs4fZz46OfZK3TYTOD7fJOm8iLgvIu5Xqly9uZCvqUoVFEn6gqTPjjc/EXFNRFzX5qvNrreLlfb1cjlPT0g6SNLnJ5Cfc/P+kqQLlW5QlJ6fJnl5sPB2FTW+VvZtX2Vfk7R3k7yUlh+M1UG9sFnZUvr+mUCdtZTzvJCvTuq5fd9uE6x/l7bNOqyvN7ue7izphIh4PCJuknR9zmtte02T9ES+Tu+pdNO2l/lqptTjv4O4o+/bq9/B8TMkLSy8v1UlBJODkJv5X6R0N6aM5S9ve76ku5QO0lLSkfR1pYv4kyUtvyYknevUDfLfS0rjWZLulvTj3F3jSNurlJRWzbslzStjwRFxm1Kr3C2S7pD0QEScW0JSV0l6te11nLqjbq/U8lGmp0bEHfn1nZKeWnJ6w6LZ8dKsrPy2pH2VLhSHKF2EDomIcZ+vHR5XS/KTg6AHJK0j6VCl4HMfSUdIOljpTvJ4dXLs9W3bFHRyfDbL13mSZigFjd+0vZOkyyKijBtbbfMTEdcolYuXSTpT0saSlqtVNHvgQ0qtSAPLj+2DbS+UtJukzzWYpW/7yvbOkm6LiMtbzDZsx05VNdsPw7J/Xpm7xf7c9hat8lzSeTVDjeu5A91uberffd9mHdTXm11Pm22v45UCwfOUrm0fl3RMRDza43xJ0ttyt/KTbdeuvWXvx6+rddzR9+01pavsoyHbq0o6RdKedXeteyYinpD0wtwacprtLSOip89T295B0l0Rcantmb1cdgPbRMRttp8i6Tzb1+Y7p700Ramb0ici4iLb31DqElnKXXfbKyh1M9mnpOWvpXTCP0vSIkk/tf3eiDi2l+lExDW2vyTpXEmPSJqvdAe1LyIibJfSM2KYjOd4iYhblLqKyfbGSi1019g+Rql72Gcj4q9d5mPcx1VEzFfqRifbr1EKrm37RKW7t3tFxN87zctEjr0ytk2TdLo6PvPF/D05X1OVutfubPtwSRsqtbicMdF8dSMi9qy9tn2mpI/Z3k+pS9x5EfGD8Sw3L+Nfko4bZH4iYj9J+9neR9Ie6rxVvKf7Kt/g2Vep23nXhvHYwVJ93j+XKT1+87Dt7SWdLmmTNvnbs/Z6oudVL+u5vdxubfI1kG3W6/p6RDyg1LJau17PkfQW2z9Q6t58WET8sQf5OlPSvIh43PbHlG40v67F8ia8H8uIO3qxvfrdcnybxrYCrJ8/G1n5gDhF0nERcWrZ6UXqFvxrpS4Nvba1pJ1sL1DqZvI62z0NvGpya5Ui4i5Jpyl1hei1WyXdWrg7drJSsFyW7ZTumnUcDHTp9ZJuioi7I2Kx0vOeryojoYj4YUS8JCJeI+l+ped5yvT33A1L+f9dJac3DFodL52UlbVW2k8qdUfaW112mcs6Oa6W5Cd3UVpD0r21ibmL0/5KXao+n/Pyg5y3rnRw7PVz29R0cnx2kq+PK7W0b6V09/tdkvaaQL5aaZuf3Kp5qaRVJW0UEe+U9HaPY3A327Ml7SBpt4hodPOgr/nJjpP0tvHkRb3ZVxsp3XS6PF9X15d0me2nDSg/aK3Zfhj4/omIB2vdYiPiHElTnQZ3K/286qCeO5Dt1i5fg9xmOc1Falxfb3Y97WR7fVbp+rarpN9J+oCkAzrNU6t8RcS9EfF4fnukpJfU57dFvsa7HzuJO/q+vfodHP9J0ia2n5VbTd4taWTveuYK4Q8lXRMRh5eYznr5To9sryTpDZKu7XU6EbFPRKwfETOU9s2vIuK9vU7H9iq2V6u9Vrqr3vNRxSPiTkkLbW+aP9pW0l96nU7BriqpS3V2i6StbK+cj71tlZ6z6bncoi/bGyo983l8GekUnKFUaCn//1nJ6Q2DVsfLGZLe72Qrpa7OtW69sv1aSbdHxN+Unq15Mv+NJ4jo5Lgq7p+3K5UNxQDo/ZLOiYj7JpqfDo69fm6bYprtjs9fSHqj7bXy3eo35s9q+VpLKXg8upCvkLTSBPLVSsvrba5Y7qn0HNZKWvr82fJKLe0ds/1mpRsQO7XoutaX/DiPmp7trMbXyr7sq4i4MiKeEhEz8nX1VqVBhO4cRH7QVrOyZeD7x/bTcvks2y9Xqr/fq5LPqw7ruX3fbp3kaxDbrMP6erPr6RmS3u00OvOzlFq5Ly4sexNJ60ca06Pb7dU2X7UbwNlOWloPKG0/dhh39H17dTWaWC/+lJ4h+6vSqGL7lZjOPKUufouVLkgfLiGNbfKGvkKp+998SduXkM7zJf05p3OVpM/1YT/NVEmjVSuNEHh5/ru65OPghZIuydvudElrlZTOKkqF7hol75cDlQq0qyQdI2laSen8VulGwuWStu3xspc5N5WeH/k/SX9TGqF47TK346D/Gh0vknaXtHt+baVnaG+QdKUKIxLnaefVtpHSIFaX5WN8614dV0oDkeyUp68o6adKA15cLOnZhe+urHQXemp+/+qc50slbdqLY6+f26ab41PSS5VHN83vP5S30fWSPli33K8pj+Kdt+e5SuXfJ8aRn7fk149L+rukX+R5pyvdpKh9t+n1VqnyN7uw3ebl7dlypPgm+ble6dmv+fnvu/3IT5O8nKJ0HF+h1E3wGYPcV3XTFyiPVt2P/PDX0fHSadlS6v7pIG975GVervSM56sK3+35eV74fsN67qC3W4f56vs2U5P6ujq/nu6X83SdpO3qln2SpE3y66dI+kNev7f1KF+HFrbXryU9t1/Hf/7uTOW4Y9Dby/lLAAAAAABUVr+7VQMAAAAAMHQIjgEAAAAAlUdwDAAAAACoPIJjAAAAAEDlERwDAAAAACqP4BgAAAAAUHkExwAAAACAyvv/NyiQ+k/A+yQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1224x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 3, figsize=(17, 4));\n",
    "axs[0].hist(sessions_per_user, bins=np.arange(0, 11)-0.5, range=(0, 11), density=True, rwidth=0.5);\n",
    "axs[0].set_title('Sessions per user distribution');\n",
    "axs[0].get_yaxis().set_visible(False);\n",
    "axs[0].set_xticks([x for x in range(0, 11)]);\n",
    "axs[0].set_xticklabels([str(x) for x in range(0, 11)]);\n",
    "\n",
    "\n",
    "axs[1].hist(baseline_conversion_rates, bins=100, density=True);\n",
    "axs[1].set_title('Baseline conversion rates');\n",
    "axs[1].xaxis.set_major_formatter(ticker.PercentFormatter(xmax=1));\n",
    "axs[1].get_yaxis().set_visible(False);\n",
    "\n",
    "\n",
    "axs[2].hist(conversion_uplifts, bins=100, density=True);\n",
    "axs[2].set_title('Converison rate uplifts (absolute)');\n",
    "axs[2].xaxis.set_major_formatter(ticker.PercentFormatter(xmax=1));\n",
    "axs[2].get_yaxis().set_visible(False);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T21:09:08.936691Z",
     "iopub.status.busy": "2021-01-31T21:09:08.936436Z",
     "iopub.status.idle": "2021-01-31T21:09:09.590692Z",
     "shell.execute_reply": "2021-01-31T21:09:09.590158Z",
     "shell.execute_reply.started": "2021-01-31T21:09:08.936650Z"
    }
   },
   "outputs": [],
   "source": [
    "data = {\n",
    "    'user': [],\n",
    "    'session_id': [],\n",
    "    'assignment': [],\n",
    "    'session_converted': []\n",
    "}\n",
    "\n",
    "# Simulate all sessions for each user\n",
    "for user_id, num_sessions in enumerate(sessions_per_user):\n",
    "    for session_id in range(1, num_sessions+1):\n",
    "        # randomly assign session to control (0) or test (1)\n",
    "        assignment = np.random.randint(0, 2)\n",
    "        \n",
    "        # if assigned to test, give them a conversion boost\n",
    "        new_conversion_rate = baseline_conversion_rates[user_id] + assignment*conversion_uplifts[user_id]\n",
    "        \n",
    "        # see if the session converted\n",
    "        session_converted = np.random.choice([0, 1], p=[1-new_conversion_rate, new_conversion_rate])\n",
    "        \n",
    "        # record the results\n",
    "        data['user'].append(user_id)\n",
    "        data['session_id'].append(f\"{user_id}-{session_id}\")\n",
    "        data['assignment'].append(assignment)\n",
    "        data['session_converted'].append(session_converted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T21:09:09.592177Z",
     "iopub.status.busy": "2021-01-31T21:09:09.592012Z",
     "iopub.status.idle": "2021-01-31T21:09:09.620540Z",
     "shell.execute_reply": "2021-01-31T21:09:09.619932Z",
     "shell.execute_reply.started": "2021-01-31T21:09:09.592156Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average sessions per user: 2.020\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>user</th>\n",
       "      <th>session_id</th>\n",
       "      <th>assignment</th>\n",
       "      <th>session_converted</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0-1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>0-2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>0-3</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>1-1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2</td>\n",
       "      <td>2-1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20191</th>\n",
       "      <td>9998</td>\n",
       "      <td>9998-2</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20192</th>\n",
       "      <td>9998</td>\n",
       "      <td>9998-3</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20193</th>\n",
       "      <td>9999</td>\n",
       "      <td>9999-1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20194</th>\n",
       "      <td>9999</td>\n",
       "      <td>9999-2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20195</th>\n",
       "      <td>9999</td>\n",
       "      <td>9999-3</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>20196 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       user session_id  assignment  session_converted\n",
       "0         0        0-1           1                  0\n",
       "1         0        0-2           0                  0\n",
       "2         0        0-3           1                  0\n",
       "3         1        1-1           1                  0\n",
       "4         2        2-1           0                  1\n",
       "...     ...        ...         ...                ...\n",
       "20191  9998     9998-2           1                  0\n",
       "20192  9998     9998-3           0                  0\n",
       "20193  9999     9999-1           1                  0\n",
       "20194  9999     9999-2           0                  0\n",
       "20195  9999     9999-3           1                  0\n",
       "\n",
       "[20196 rows x 4 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame(data)\n",
    "print (f\"Average sessions per user: {df.shape[0] / len(df.user.unique()):0.3f}\")\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T21:09:09.621816Z",
     "iopub.status.busy": "2021-01-31T21:09:09.621651Z",
     "iopub.status.idle": "2021-01-31T21:09:09.641646Z",
     "shell.execute_reply": "2021-01-31T21:09:09.641065Z",
     "shell.execute_reply.started": "2021-01-31T21:09:09.621793Z"
    }
   },
   "outputs": [],
   "source": [
    "num_samples = 50000\n",
    "\n",
    "control_converted = df[df.assignment==0].session_converted.sum()\n",
    "control_total = df[df.assignment==0].session_converted.count()\n",
    "control_samples = np.random.beta(control_converted, control_total - control_converted, size=num_samples)\n",
    "\n",
    "\n",
    "test_converted = df[df.assignment==1].session_converted.sum()\n",
    "test_total = df[df.assignment==1].session_converted.count()\n",
    "test_samples = np.random.beta(test_converted, test_total - test_converted, size=num_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T21:09:09.642945Z",
     "iopub.status.busy": "2021-01-31T21:09:09.642617Z",
     "iopub.status.idle": "2021-01-31T21:09:09.646943Z",
     "shell.execute_reply": "2021-01-31T21:09:09.646057Z",
     "shell.execute_reply.started": "2021-01-31T21:09:09.642913Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test converts higher than control: 100% of the time\n"
     ]
    }
   ],
   "source": [
    "test_gt_control = (test_samples > control_samples).mean()\n",
    "print(f\"Test converts higher than control: {test_gt_control:0.0%} of the time\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T21:09:09.648355Z",
     "iopub.status.busy": "2021-01-31T21:09:09.648052Z",
     "iopub.status.idle": "2021-01-31T21:09:10.510653Z",
     "shell.execute_reply": "2021-01-31T21:09:10.509958Z",
     "shell.execute_reply.started": "2021-01-31T21:09:09.648326Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHQAAAF1CAYAAACJaL20AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABW/0lEQVR4nO3deXxU1f3/8fchggmgCARQDC2CbCaEsEMEZC2CAi644BpxKbW18lVbF4ryrdqiP1TUSqmKX9BaBEVEUKyipMguS0CQRVQUSkSIgLIneH5/3Js4hJnJzr138no+HvOYycydO5/zuZOZM5977rnGWisAAAAAAAAERxWvAwAAAAAAAEDJUNABAAAAAAAIGAo6AAAAAAAAAUNBBwAAAAAAIGAo6AAAAAAAAAQMBR0AAAAAAICAoaADeMgY090Ys8nrOPzKGLPfGNPE6zgAAKgsjDETjTGjQ/7+jTFmp/udXNcYc74x5nP370s8DLXUjDHrjTE9vY6jKMaYa40x73sdhx8ZY37hvgfjvI4F8JKx1nodA+Arxphukh6XlCzpmKQNkkZaaz/xNLASMsZUk/SApGslNZS0S9JHkv5srd3qYWieM8ZsldRAzvbdL+k9Sb+z1u4vxnMzJN1ire1WkTECAFDeQr7/8uR8B34m6WVJz1trfwqzfFVJP0jqYq1d4973oaS3rbVPn6y4/cQYM0TS/0pqIumopLWSbrbWfuVpYCVkjDld0p8lXSapjqSdkmZLesRau9vL2LxmjLGSDkqykvZJmibpD9baY8V47hhJ51prr6vQIAEXI3SAEO6X2xxJz8r5cjtbzpf2ES/jKqU3JA2WdI2kWpLaSFopqY+XQYUyxpzi4csPstbWlJQmqa2k+z2MBQCAk2WQtfY0Sb+UNFbSvZImRVi2gaR4SetD7vtlob+LzePv/TIzxpwrpwB2t5y+1TmSnpNTHAsMd6ffh3J2Xl4o6XRJXSXlSOrkYWjH8fj90sbtJ14g6SpJwz2MBYiIgg5wvOaSZK2daq09Zq09ZK1931q7Nn8BY8xwY8wGY8weY8y/jTG/dO83xpinjDHfGWN+MMZ8aoxJcR8baIz5zBjzozHmv8aYe9z7expjtoesu5UxJtMYs9cdDjw45LHJxpjnjDHvuOtZZoxpGq4Rxpi+kvpJGmKt/cRam2et3Wetfc5aO8ldpqEx5m1jzPfGmC3GmFtDnj/GGDPdGPOy+1rrjTEd3MfuNca8Uej1njbGPOPermWMmWSMyXbb+kj+cFhjTIYxZpGbpxxJY4wx5xpj/mOM2WeM2W2MmRayXut2nvLX+7IxZpcx5mtjzJ+MMVVC1rvQGDPO3S5fGWMGFGeDW2u/lfRvOYWd/Ne9zxjzhdv2z4wxl+ZvH0kTJXV1h/nude8/1X3tb9xh6RONMQnFeX0AALzg9gvelvNj9caQPstk97u7uaT8w8L3GmM+MsZ8IWdkymz3e/DUUnzvR/zOzO8XGWPudvtT2caYm/JjNsYkGGOecPsB+9zv/vzndjHGLHb7UGtMlEOqjDFb3b5S1D5PGGmSvrLWfmgdP1prZ1hrv3HXVSWkD5HjrreO+1i8Meaf7v17jTGfGGMahOTpS/f1vzLGXBty/8KQuNPd5+1zr9NDHss0xjzs5vtHY8z7xpjECO24QdIvJF1qrf3MWvuTtfY7a+3D1tp33fWVqk9qjPm7MWZcoXzPMsbc5d5uaIyZ4fbnvjLG/D5kuTHGmDfcPP0gKcMY08kYs8I4feudxpgn3WUbG6efeErIekvcry2KtXaLpEU6vp/4tDFmmxvTSmNMd/f+C+WMjr/K/f/IH9UW8X8EKCsKOsDxNks6ZoyZYowZYIypHfqgcYbZPiBneGo9SR9Lmuo+/CtJPeQUhWpJulLOng7J2fP1a3ePWIqcQ5+OY5xhzbMlvS+pvqQ7JL1qjGkRstjVckYM1Za0RdKjEdrRV9Jya+22KG19TdJ2OYdjDZX0F2NM75DHB7vLnCHpbUl/C3neQGPMaW7ccW5b/+U+PlnOUO5z5Yx8+ZWkW0LW21nSl3L2+j0q6WG3zbUlJckZHRXOs3Ly2kTO3pIbJN0U8nhnOR3PRDmHzE0yxpgo7Zcbf5KkAXLyme8LSd3d1/tfSf80xpxlrd0gaYSkJdbamtbaM9zlx8rZ7mluu8+W9GBRrw0AgNestcvl9Ae6F7p/s5wRHJJ0hrW2t7W2qaRv5I5ytdYeUcm/94v6zjxTzvfv2ZJulvRcSH9snKT2ktLljKT+o6SfjDFnS3pH0iPu/fdImmGMqVfMNETq8xS2SlJL4xSoehljahZ6/A5Jl8jppzSUtEfOCB5JutFtVyNJdeX0Jw4ZY2pIekbSALefmC4pq/ALu4Whd9xl60p6UtI7xpi6IYtdI6dvVF9SNTl5CKevpPciHWpexj7pVDkFDeOuq7ac98RrxtkRN1vSGjnbt4+kkcaY/iHrHSJnlPkZkl6V9LSkp621p0tqKml6hDaVtl8blTGmpZz/jdB+4idy3r915PR/XzfGxFtr35P0F0nT3P+PNu7ykxX9fwQoNQo6QAhr7Q+Susk5ZvYFSbvcan8Dd5ERkv5qrd1grc2T86GdZpxROrmSTpPUUs78VBustdnu83IlnWeMOd1au8dauyrMy3eRVFPSWGvtUWvtR3IO/xoWssxMa+1y97VfVcjegkLqSsqO8JiMMY0knS/pXmvtYWttlqQX5RRJ8i201r7rHi/8ipxDtmSt/VpOh+ZSd7nekg5aa5e6eRooZ86hA9ba7yQ9JedLP98Oa+2z7qihQ25ufimpoRvLQhXiFo2ulnS/uzdsq6QnJF0fstjX1toX3HinSDpLTucxkreMMT9K2ibpO0kP5T9grX3dWrvD3WM1TdLnijAE2e2w3Cbpf6y131trf5Tzvrg63PIAAPjQDjk/TkukpN/7kg6r6O/MXDnz/eW6o0X2S2rhFgOGS7rTWvtfdyT1YreodJ2kd91+y0/W2g8krXBjK46wfZ7CrLVfSuoppxgxXdJud7RKfmFnhKRR1trtblxjJA11R5HkyumfnevGvtLtd0rST5JSjDEJ1tpsa224Q9oukvS5tfYVtw81VdJGSYNClvk/a+1mt381XaXsJ6psfdKP5fSj8wuEQ+XsCNshqaOketbaP7vr/VJOfzt0+y+x1r7lbsf8fuK5xphEa+1+a+3SwsGWpV8bxSpjzAE5c2lmSpqQ/4C19p/W2hx3Ozwh6VRJLcKtpJj/I0CpUdABCnELMRnW2iQ5o2kaShrvPvxLSU+7w0/3SvpekpF0tvtl9zc5e2K+M8Y8b5w5eSTpcjkf5l8b5/CirmFeuqGkbfb4SQm/ltNpyPdtyO2Dcr5sw8mRU9CIpKGk/I5UcV8r3vx8LPO/9POX+jX6eXTOLyVVlZQdkqN/yNm7k6/wqKE/ysnhcncIbLhjlBPd9X5dnHittQfdm5HyI0mXuHvCesopwhUMSzbG3GCMyQppQ0ro44XUk1Rd0sqQ5d9z7wcAIAjOltOnKamSfu8X5zszxy0S5Mvv7yTKmc/niwhxXJG/Tne93RS9LxQqWp/nONbapdbaK6219eQULXpIGhUSx8yQGDbImV+ngZwiwr/ljFTZYYx53BhT1Vp7QM5hbyPk5PEdd1RIYQ11fD9Iqth+Yqn6pNZaK2ckTGg/8VX39i8lNSy0nR7Q8TvgCvcTb5YzomujcQ4zuzhCvGXp14bTzm3TVXJGmdXIf8AYc49xpl/Y57ahliL3E4vzPwKUGgUdIApr7UY5wyRT3Lu2yTl06oyQS4K1drG7/DPW2vaSzpPz5fMH9/5PrLVD5Hx4v6Xww0V3SGrk7oHK9wtJ/y1F6PMkdXIPJwpnh6Q6+YdNleK1XpfU013/pfq5oLNNzgTSiSH5Od1amxzy3ONOrWet/dZae6u1tqGkX0uaYNx5c0Ls1s8jeUoTb0TW2v/I2cbjJMkdbfWCpN9Jqmudw6rWySk6nRC/G9shSckhba5lnYn0AADwNWNMRzk/fE8YIVsMJf3eL8t35m45I3zCzR+4TdIrhfpnNay1Y0vRpmKzzhlQ39Tx/cQBheKId0cU5Vpr/9dae56cw6ouljuCxFr7b2ttPzlFlo1y+iGF7dDx/SCpbP3E/u7hXuGUtU86Vc7IpF/KKYbMcO/fJmcOotD8nGatDR1JVbif+Lm1dpicPvRjkt4IE3dZ+7VhWcd0SUvkHhbozpfzRznTDdR2+4n7FLmfWJz/EaDUKOgAIYwxLY0zEV+S+3cjOXsY8od3TpR0vzEm2X28ljHmCvd2R2NMZ/e44wNyOh0/GWOqGWOuNcbUstbmyjn95wmnBpW0TM4egz8aY6oaZzK/QXL2cpSItXaepA/k7CVqb4w5xRhzmjFmhDFmuHXm1lks6a/GmaQvVc4ekH8Wc/275Aw//T85X8wb3Puz5Rxv/YQx5nTjTA7Y1BhzQaR1GWOuCCk87ZHzRXhcftzhsdMlPeq245eS7ipuvMUwXlI/Y0wbOXtgrJzTvMs4kzGmhCy7U1KScc4QIXfv1QuSnjLG1Hefc3ah48EBAPAV93v6Yjn9jH9aaz8t6TpK+r1flu9M97kvSXrSOBPgxhljuhpjTpXTHxhkjOnv3h9vnAmWI+3YKhVjTDdjzK0hsbeUMzdLaD/xUfPzCTPqGWf+RRlnzp3WxjmM/Ac5O6p+MsY0MMYMcYsUR+QcYhaun/iupObGmGvcft1VcnYgzilFU16RU2iY4fZ9qxhj6hpjHjDGDFQZ+6TW2tVyCnAvSvq3tXav+9ByST8a5wQbCe62SnGLimEZY64zxtRzt3/+egr3E8vUry2GsZJuNcacKWd6hTw5/cRTjDEPyjlLWL6dkhrnF8NK0zcGSoKCDnC8H+XsSVhmnONml8oZnXG3JFlrZ8rZO/CacWbfXydnQl3J+TB/QU5R4ms5w1n/n/vY9ZK2us8ZIenawi9srT0q58tygJwvwQmSbnBHCZXGUDlf/tPk7DlYJ6mDnL0yklOoaixnr8ZMSQ+5haDi+pecSfX+Vej+G+RMxPeZnFy8oejDejvKyfd+OZPU3ekeU13YHXIKZV/K2Yv4LzkduzJzC1QvS3rQWvuZnPl5lsj5Um4t5+wG+T6Sc7rWb40xu9377pUzWd5SdxvPU4RjqQEA8Nhs8/MccqPkTK57U/SnRFXS7/2yfGfeI+lTOZPSfi+nT1bF/UGff+KKXXLa9geV/2+dvXIKOJ+6/Zb35PShHncff1pOX+Z9N8dL5fQrJWey5zfkFHM2SPqPnMJKFTk7qXa4bbpA0m8Kv7C1NkfOqJ675fQx/yjpYmvt7sLLFsU68/v0lTMa6AM3puVyDhtaVk590hP6ie4Ouovlni1MPxd9akVZz4WS1rv5flrS1daZW6ewsvZrI3KLnQvkvKf+LWe7b5bT3z+s4w8Te929zjHG5M+ZWdL/EaDYjHOYIwAAAAAAAIKCEToAAAAAAAABQ0EHAAAAAAAgYCjoAAAAAAAABAwFHQAAAAAAgIChoAMAAAAAABAwp5THShITE23jxo3LY1UAAMCHVq5cudtaW8/rOHA8+mColH76ybmuwr7pSH5yc1SFHIX1k3XzY8gP/C9aH6xcCjqNGzfWihUrymNVAADAh4wxX3sdA05EHwwAgNgWrQ9GSRIAAABAcEyY4FwQ0YQJEzSBHEU04ZMJmvAJ+UHwUdABAAAAEBzTpzsXRDR9+nRNJ0cRTV8/XdPXkx8EHwUdAAAAAACAgCmXOXQAAPBCbm6utm/frsOHD3sdSsyIj49XUlKSqlat6nUoAACcNPQp4LXS9MEo6AAAAmv79u067bTT1LhxYxljvA4n8Ky1ysnJ0fbt23XOOed4HQ4AACcNfQp4qbR9MA65AgAE1uHDh1W3bl06XuXEGKO6deuydxIAUOnQp4CXStsHY4QOACDQ6HiVL/IJwPcyM72OwPcyyVFUmRmZYe/nOxBeKs37jxE6AACUwbfffqurr75aTZs2Vfv27TVw4EBt3ry5xOsZP368Dh48WOLn1axZs8TPAQAA/jN8+HDVr19fKSkpx93//fffq1+/fmrWrJn69eunPXv2SHIKd4sXLy5YLiMjQ2+88UaRrxMXF6e0tDSlpKToiiuuKFX/o6wKxz5x4kS9/PLLkorfjlDjx48veL4kPfvss2rZsqWSk5P1xz/+seD+tWvXqmvXrkpOTlbr1q3Djohp3Lixdu/eXdImHefo0aPq0aOH8vLyyrSeojBCBwAQMwYNKt/1zZ4d/XFrrS699FLdeOONeu211yRJa9as0c6dO9W8efMSvdb48eN13XXXqXr16ic8duzYMcXFxZVofQAQs8aNc67vucfbOHxsnJuje8hRWOMWu/lJ91d+MjIy9Lvf/U433HDDcfePHTtWffr00X333aexY8dq7Nixeuyxx5SZmamaNWsqPT29RK+TkJCgrKwsSdK1116riRMn6q677iryeXl5eTrllPIpIRSOfcSIEaVeV15enl566SWtWrVKkjR//nzNmjVLa9as0amnnqrvvvuuYLnrrrtOr7zyitq0aaOcnJwKOwlEtWrV1KdPH02bNk3XXntthbyGxAgdAABKbf78+apatepxnZA2bdqoW7du+sMf/qCUlBS1bt1a06ZNk+R0Xnr27KmhQ4eqZcuWuvbaa2Wt1TPPPKMdO3aoV69e6tWrlyRn5M3dd9+tNm3aaMmSJXryySeVkpKilJQUjR8/3ovmAoA/zJnjXBDRnDlzNIccRTRn8xzN2ey//PTo0UN16tQ54f5Zs2bpxhtvlCTdeOONeuutt7R161ZNnDhRTz31lNLS0vTxxx9LkhYsWKD09HQ1adKkWKNcunfvri1btujAgQMaPny4OnXqpLZt22rWrFmSpMmTJ2vw4MHq3bu3+vTpo/379+umm25S69atlZqaqhkzZkiS3n//fXXt2lXt2rXTFVdcof3790tyRrs89NBDateunVq3bq2NGzeGjX3MmDEFhchQK1eu1AUXXKD27durf//+ys7OPmGZjz76SO3atSsoNv3973/Xfffdp1NPPVWSVL9+/YIYU1NT1aZNG0lS3bp1I+4we/zxx9W6dWt16tRJW7ZskeQU3EaMGKEOHTqoefPmBf9j69evV6dOnZSWlqbU1FR9/vnnkqRLLrlEr776apHboCwYoQMAQCmtW7dO7du3P+H+N998U1lZWVqzZo12796tjh07qkePHpKk1atXa/369WrYsKHOP/98LVq0SL///e/15JNPav78+UpMTJQkHThwQJ07d9YTTzyhlStX6v/+7/+0bNkyWWvVuXNnXXDBBWrbtu1JbS8AAJVGz54n3nflldLtt0sHD0oDB574eEaGc9m9Wxo69PjHyjCv0c6dO3XWWWdJks4880zt3LlTjRs31ogRI1SzZs2CkViTJk1Sdna2Fi5cqI0bN2rw4MEaWjiOEHl5eZo7d64uvPBCPfroo+rdu7deeukl7d27V506dVLfvn0lSatWrdLatWtVp04d3XvvvapVq5Y+/fRTSdKePXu0e/duPfLII5o3b55q1Kihxx57TE8++aQefPBBSVJiYqJWrVqlCRMmaNy4cXrxxRdPiP3DDz88Ib7c3FzdcccdmjVrlurVq6dp06Zp1KhReumll45bbtGiRcf1xzZv3qyPP/5Yo0aNUnx8vMaNG6eOHTtq8+bNMsaof//+2rVrl66++urjDscKld/Gl19+WSNHjiwo3mzdulXLly/XF198oV69emnLli2aOHGi7rzzTl177bU6evSojh07JklKSUnRJ598UsTWLRsKOgAAlLOFCxdq2LBhiouLU4MGDXTBBRfok08+0emnn65OnTopKSlJkpSWlqatW7eqW7duJ6wjLi5Ol19+ecH6Lr30UtWoUUOSdNlll+njjz+moAMAQCVjjIk6ee4ll1yiKlWq6LzzztPOnTvDLnPo0CGlpaVJckbo3HzzzUpPT9fbb79dMErm8OHD+uabbyRJ/fr1Kxg5NG/evILDzCWpdu3amjNnjj777DOdf/75kpz5Y7p27VqwzGWXXSZJat++vd58881it3XTpk1at26d+vXrJ8k5BD2/sBUqOztbrVq1Kvg7Ly9P33//vZYuXapPPvlEV155pb788kvl5eVp4cKF+uSTT1S9enX16dNH7du3V58+fU5Y57Bhwwqu/+d//qfg/iuvvFJVqlRRs2bN1KRJE23cuFFdu3bVo48+qu3bt+uyyy5Ts2bNJDl9uWrVqunHH3/UaaedVux2lwQFHQAASik5ObnEk/blD/+VnC/6SJPlxcfHM28OAABeiTaipnr16I8nJpbr2dgaNGig7OxsnXXWWcrOzi44hCic0H6GtTbsMqFz6IQuO2PGDLVo0eK4+5ctW1awQykSa6369eunqVOnRo0pWr8n0nqTk5O1ZMmSqMslJCQcN7lxUlKSLrvsMhlj1KlTJ1WpUkW7d+9WUlKSevToUTAaeuDAgVq1alXYgk5o0SzS7fy/r7nmGnXu3FnvvPOOBg4cqH/84x/q3bu3JOnIkSOKj48vdptLijl0UCkMGvTzBQDKS+/evXXkyBE9//zzBfetXbtWZ5xxhqZNm6Zjx45p165dWrBggTp16hR1Xaeddpp+/PHHsI91795db731lg4ePKgDBw5o5syZ6t69e7m2BQACIyHBuUDS8f3c/EtCQoISyFFECVUTlFA1OPkZPHiwpkyZIkmaMmWKhgwZIil636Gk+vfvr2effbagCLR69eqwy/Xr10/PPfdcwd979uxRly5dtGjRooK5Zg4cOFDkGT+LE3uLFi20a9eugoJObm6u1q9ff8JyrVq1KnhtyRmlNH/+fEnO4VdHjx5VYmKi+vfvr08//VQHDx5UXl6e/vOf/+i8884L+9r58x9OmzbtuNFGr7/+un766Sd98cUX+vLLL9WiRQt9+eWXatKkiX7/+99ryJAhWrt2rSQpJydHiYmJFTbxskRBBwCAUjPGaObMmZo3b56aNm2q5ORk3X///brmmmsKJt3r3bu3Hn/8cZ155plR13XbbbfpwgsvLJgUOVS7du2UkZGhTp06qXPnzrrllls43ApA5TV3rnNBRHPnztVcchTR3Gvnau61/svPsGHD1LVrV23atElJSUmaNGmSJOm+++7TBx98oGbNmmnevHm67777JEmDBg3SzJkzj5sUubRGjx6t3NxcpaamKjk5WaNHjw673J/+9Cft2bNHKSkpatOmjebPn6969epp8uTJGjZsmFJTU9W1a1dt3Lgx6usVJ/Zq1arpjTfe0L333qs2bdooLS3tuFOd5xswYIAWLFhQ8Pfw4cP15ZdfKiUlRVdffbWmTJkiY4xq166tu+66Sx07dlRaWpratWuniy66KOxr79mzR6mpqXr66af11FNPFdz/i1/8Qp06ddKAAQM0ceJExcfHa/r06UpJSVFaWprWrVtXcJay+fPnR1x/eTGRhmGVRIcOHeyKFSvKIRygYoSOzCnqNMQAgmPDhg3HHTON8hEur8aYldbaDh6FhAjogwEINwKd/m7J0acItksvvVSPP/54wfw1FSEjI0MXX3xx1ImmQ1122WUaO3asmjdvXuzXKGkfjBE6AAAAAILj4YedCyJ6+OGH9TA5iujh/zysh/9DfmLJ2LFjw57S3CtHjx7VJZdcUqJiTmkwKTIAAACA4Mg/vXGEQ0Lw8ymg8w+bYRTP8T78ys3PBbyHYkWLFi1OmNC5vE2ePLnYy1arVq3g0KuKxAgdAAAAAACAgKGgAwAAAAAAEDAUdAAAAAAAAAKGOXQQs8IdKwwAAICAq1vX6wh8ry45iqpudfKD2MAIHQAASiknJ0dpaWlKS0vTmWeeqbPPPrvg76NHj5ZqnZmZmVq8eHE5R+rd6wBAuZsxw7kgohkzZmgGOYpoxpUzNONK/+XHGKPrrruu4O+8vDzVq1dPF198cYW/9siRI7VgwQJJ0t/+9jede+65MsZo9+7dBcu8+uqrSk1NVevWrZWenq41a9ZIkg4fPqxOnTqpTZs2Sk5O1kMPPRT2NXr27KkVK1aUOda+fftqz549ZV5PLGCEDgAgdpT30LwiTgFSt25dZWVlSZLGjBmjmjVr6p577il4PC8vT6ecUrKv2szMTNWsWVPp6eklDtePrwMAAIqnRo0aWrdunQ4dOqSEhAR98MEHOvvssyv8dXNycrR06VKNHz9eknT++efr4osvVs+ePY9b7pxzztF//vMf1a5dW3PnztVtt92mZcuW6dRTT9VHH32kmjVrKjc3V926ddOAAQPUpUuXCon3+uuv14QJEzRq1KgKWX+QMEIHAIBylJGRoREjRqhz58764x//qC+++EIXXnih2rdvr+7du2vjxo2SpNmzZ6tz585q27at+vbtq507d2rr1q2aOHGinnrqKaWlpenjjz9WRkaGfvOb36hLly5q0qSJMjMzNXz4cLVq1UoZGRkFr/v++++ra9euateuna644grt379fktS4cWM99NBDateunVq3bq2NGzeGfR0ACIz773cuiOjcc+/Xuefer0GDmIYgnPvn3a/75/nzPTRw4EC98847kqSpU6dq2LBhBY8dOHBAw4cPV6dOndS2bVvNmjVLkrR161Z1795d7dq1U7t27QpG4GZmZqpnz54aOnSoWrZsqWuvvVbW2hNec8aMGbrwwgsL/m7btq0aN258wnLp6emqXbu2JKlLly7avn27JGdkUc2aNSVJubm5ys3NlTEmbPteeeUVpaWlKSUlRcuXL5fk7BS7/vrr1bVrVzVr1kwvvPCCJCk7O1s9evQoWD6/vzJ48GBNnTq1mBmNbYzQQaUW+gVXxI54ACi27du3a/HixYqLi1OfPn00ceJENWvWTMuWLdPtt9+ujz76SN26ddPSpUtljNGLL76oxx9/XE888YRGjBhx3EifSZMmac+ePVqyZInefvttDR48WIsWLdKLL76ojh07KisrS0lJSXrkkUc0b9481ahRQ4899piefPJJPfjgg5KkxMRErVq1ShMmTNC4ceP04osvnvA6ABAYS5Z4HYHv7dlTdI4KF3oqU194yfai81N4dIokXXnllbr99tt18OBBDRw48ITHMzIylJGRod27d2vo0KHHPZaZmVms2K6++mr9+c9/1sUXX6y1a9dq+PDhBYWMRx99VL1799ZLL72kvXv3qlOnTurbt6/q16+vDz74QPHx8fr88881bNiwgkObVq9erfXr16thw4Y6//zztWjRInXr1u2411y0aNEJ8RZl0qRJGjBgQMHfx44dU/v27bVlyxb99re/VefOncM+7+DBg8rKytKCBQs0fPhwrVu3TpK0du1aLV26VAcOHFDbtm110UUXaerUqerfv79GjRqlY8eO6eDBg5Kk2rVr68iRI8rJyan080VR0AEAoJxdccUViouL0/79+7V48WJdccUVBY8dOXJEklP0ueqqq5Sdna2jR4/qnHPOibi+QYMGyRij1q1bq0GDBmrdurUkKTk5WVu3btX27dv12Wef6fzzz5ckHT16VF27di14/mWXXSZJat++vd58881yby8AACgfqamp2rp1q6ZOnXpC0ej999/X22+/rXHjxkly5q755ptv1LBhQ/3ud79TVlaW4uLitHnz5oLndOrUSUlJSZKktLQ0bd269YSCTnZ2turVq1fsGOfPn69JkyZp4cKFBffFxcUpKytLe/fu1aWXXqp169YpJSXlhOfmjzjq0aOHfvjhB+3du1eSNGTIECUkJCghIUG9evXS8uXL1bFjRw0fPly5ubm65JJLlJaWVrCe+vXra8eOHRR0vA4AAIBYU6NGDUnSTz/9pDPOOKNgnp1Qd9xxh+666y4NHjxYmZmZGjNmTMT1nXrqqZKkKlWqFNzO/zsvL09xcXHq169fxOHH+c+Ji4tTXl5eKVsFAEDlEW1ETfXq1aM+npiYWOwROeEMHjxY99xzjzIzM5WTk1Nwv7VWM2bMUIsWLY5bfsyYMWrQoIHWrFmjn376SfHx8QWPhfYbIvUDEhISdPjw4WLFtnbtWt1yyy2aO3du2GLKGWecoV69eum9994LW9ApfChW/t/h7u/Ro4cWLFigd955RxkZGbrrrrt0ww03SHKKWQkJCcWKOZYxhw4AABXk9NNP1znnnKPXX39dktMRyz8jxL59+womOpwyZUrBc0477TT9+OOPJXqdLl26aNGiRdqyZYsk5xj70L1z4ZTmdQAAQMUbPny4HnrooYIRufn69++vZ599tmAenNWrV0ty+hRnnXWWqlSpoldeeUXHjh0r0eu1atWqoA8RzTfffKPLLrtMr7zyipo3b15w/65duwpG2hw6dEgffPCBWrZsGXYd06ZNkyQtXLhQtWrVUq1atSRJs2bN0uHDh5WTk6PMzEx17NhRX3/9tRo0aKBbb71Vt9xyi1atWiXJ6U99++23Yef5qWwo6AAAUIFeffVVTZo0qeBUnvkTGI4ZM0ZXXHGF2rdvr8TExILlBw0apJkzZ5ZosuJ69epp8uTJGjZsmFJTU9W1a9eCyZcjKc3rAIAvJCU5F0SUkJCkhARyFEnS6UlKOt2/+UlKStLvf//7E+4fPXq0cnNzlZqaquTkZI0ePVqSdPvtt2vKlClq06aNNm7cWDBSuLguuuii40YUPfPMM0pKStL27duVmpqqW265RZL05z//WTk5Obr99tuVlpamDh06SHIO2erVq5dSU1PVsWNH9evXL+Kp1uPj49W2bVuNGDFCkyZNKrg/NTVVvXr1UpcuXTR69Gg1bNhQmZmZatOmjdq2batp06bpzjvvlCStXLlSXbp0KfGZRGORCTfLdUl16NDBlsf55IHyFGlG/9AJ35gUGeWKN9RJt2HDBrVq1crrMGJOuLwaY1Zaazt4FBIioA8GoLzOYlXZuy6VvU/RrVs3zZkzR2ecccZJf+0xY8aU6EQNd955pwYPHqw+ffpUcGQnX0n7YIzQAQAAAACgEnviiSf0zTffeB1GsaSkpMRkMac0GKMEAAAAIDhGjnSux4/3MgpPFHc0zvr1IyVJycnjKyyWIBv53khJ0vgLx3sah59EOs34yRDtxBDh3HrrrRUTSABR0EGlU17DUgEAAOCBMGcOjFWl7bfu25dVrnHEmqxvs7wOASgXFHQQOIW/2Cr78b7wGBVCz1lrTzjVJUqvPObWAwAAQMWjoINAOBm/mZnPFgie+Ph45eTkqG7duhR1yoG1Vjk5OYqPj/c6FACIeeykBFBWFHQQUxgsAVQu+afU3LVrl9ehxIz4+HglcTpgAAAA36OgAwAIrKpVq+qcc87xOgwAwMnUvLnXEfhezZrkKJrmdf2Zn6efflovvPCCrLW69dZbNdKdAPz777/XVVddpa1bt6px48aaPn26ateurczMTFWrVk3p6emSpIyMDF188cUaOnRo1NeJi4tT69atlZeXp1atWmnKlCmqXr16RTfvOIVjnzhxoqpXr64bbrih2O0INX78eNWpU0c33HCDXn/9dY0ZM0YbNmzQ8uXL1aGDc8bv3Nxc3XLLLVq1apXy8vJ0ww036P7775ckvffee7rzzjt17Ngx3XLLLbrvvvtOeI2ePXtq3LhxBesrrb59++r1119X7dq1y7QeidOWIwYMGvTzBQAAADHu+eedCyJKTX1eqankKJLnBz2v5wf5Kz/r1q3TCy+8oOXLl2vNmjWaM2eOtmzZIkkaO3as+vTpo88//1x9+vTR2LFjJTlFkcWLF5f4tRISEpSVlaV169apWrVqmjhxYrGel5eXV+LXiqRw7CNGjNANN9xQqnXl5eXppZde0jXXXCPJOa35m2++qR49ehy33Ouvv64jR47o008/1cqVK/WPf/xDW7du1bFjx/Tb3/5Wc+fO1WeffaapU6fqs88+K33jinD99ddrwoQJ5bIuCjoAAAAAAHhow4YN6ty5s6pXr65TTjlFF1xwgd58801J0qxZs3TjjTdKkm688Ua99dZb2rp1qyZOnKinnnpKaWlp+vjjjyVJCxYsUHp6upo0aaI33nijyNft3r27tmzZogMHDmj48OHq1KmT2rZtq1mzZkmSJk+erMGDB6t3797q06eP9u/fr5tuukmtW7dWamqqZsyYIUl6//331bVrV7Vr105XXHGF9u/fL0lq3LixHnroIbVr106tW7fWxo0bw8Y+ZswYjRs37oT4Vq5cqQsuuEDt27dX//79lZ2dfcIyH330kdq1a6dTTnEOQGrVqpVatGhxwnLGGB04cEB5eXk6dOiQqlWrptNPP13Lly/XueeeqyZNmqhatWq6+uqrC9pf2CuvvKK0tDSlpKRo+fLlkpzTrl9//fXq2rWrmjVrphdeeEGSlJ2drR49ehQsn7+NBg8erKlTpxa5bYqDQ64AoKQYDgYAgHduu825ZpRORGvXOjlilE54t8128hNtlE7PyT1PuO/K5Ct1e8fbdTD3oAa+OvCExzPSMpSRlqHdB3dr6PTjDxfKzMiMGlNKSopGjRqlnJwcJSQk6N133y04tGfnzp0666yzJElnnnmmdu7cqcaNG2vEiBGqWbOm7rnnHknSpEmTlJ2drYULF2rjxo0aPHhw1MOW8vLyNHfuXF144YV69NFH1bt3b7300kvau3evOnXqpL59+0qSVq1apbVr16pOnTq69957VatWLX366aeSpD179mj37t165JFHNG/ePNWoUUOPPfaYnnzyST344IOSpMTERK1atUoTJkzQuHHj9OKLL54Q+4cffnhCfLm5ubrjjjs0a9Ys1atXT9OmTdOoUaP00ksvHbfcokWL1L59+6j5laShQ4dq1qxZOuuss3Tw4EE99dRTqlOnjv773/+qUaNGBcslJSVp2bJlYddx8OBBZWVlacGCBRo+fLjWrVsnSVq7dq2WLl2qAwcOqG3btrrooos0depU9e/fX6NGjdKxY8d08OBBSVLt2rV15MiRghN7lAUFHQAIh1NPAADgT5s3ex2B7+3fT46i2Zzjv/y0atVK9957r371q1+pRo0aSktLU1xc3AnLGWOintnzkksuUZUqVXTeeedp586dYZc5dOiQ0tLSJDkjdG6++Walp6fr7bffLhglc/jwYX3zzTeSpH79+qlOnTqSpHnz5um1114rWFft2rU1Z84cffbZZzr//PMlSUePHlXXrl0LlrnsssskSe3bty8YdVQcmzZt0rp169SvXz9J0rFjxwoKW6Gys7PVqlWrIte3fPlyxcXFaceOHdqzZ4+6d+9eULQqrmHDhkmSevTooR9++EF79+6VJA0ZMkQJCQlKSEhQr169tHz5cnXs2FHDhw9Xbm6uLrnkkoKcS1L9+vW1Y8cOCjoAAAAAgLIJNwC5Mu/PijaipnrV6lEfT6yeWOSInHBuvvlm3XzzzZKkBx54oOCskw0aNFB2drbOOussZWdnq379+hHXceqppxbcttaGXSZ/Dp1Q1lrNmDHjhEOVli1bpho1akSN21qrfv36RTyMKD+muLi4Es3DY61VcnKylixZEnW5hIQEHT58uMj1/etf/9KFF16oqlWrqn79+jr//PO1YsUKNWrUSNu2bStYbvv27Tr77LPDrqNwMS3/73D39+jRQwsWLNA777yjjIwM3XXXXQXzBB0+fFgJCQlFxlwU5tCBbzHZMU4K3mgAAADwge+++06S9M033+jNN98smOR38ODBmjJliiRpypQpGjJkiCTptNNO048//lgur92/f389++yzBUWg1atXh12uX79+eu655wr+3rNnj7p06aJFixYVTOJ84MABbS5iJF1xYm/RooV27dpVUNDJzc3V+vXrT1iuVatWBa8dzS9+8Qt99NFHBTEuXbpULVu2VMeOHfX555/rq6++0tGjR/Xaa69p8ODBYdcxbdo0SdLChQtVq1Yt1apVS5Izz9Hhw4eVk5OjzMxMdezYUV9//bUaNGigW2+9teDsWpJTqPr222/VuHHjImMuCgUdAAAAAAA8dvnll+u8887ToEGD9Nxzz+mMM86QJN1333364IMP1KxZM82bN6/glNqDBg3SzJkzj5sUubRGjx6t3NxcpaamKjk5WaNHjw673J/+9Cft2bNHKSkpatOmjebPn6969epp8uTJGjZsmFJTU9W1a1dt3Lgx6usVJ/Zq1arpjTfe0L333qs2bdooLS0t7Fm9BgwYoAULFhT8PXPmTCUlJWnJkiW66KKL1L9/f0nSb3/7W+3fv1/Jycnq2LGjbrrpJqWmpuqUU07R3/72N/Xv31+tWrXSlVdeqeTk5LAxxcfHq23bthoxYoQmTZpUcH9qaqp69eqlLl26aPTo0WrYsKEyMzPVpk0btW3bVtOmTdOdd94pyZnouUuXLgWTOJeFiTQMqyQ6dOhgV6xYUeb1AKG8HDBRmYeXVjqhb7TQDR9tDp3ivDl5EyHGGGNWWms7eB0HjkcfDJXSyJHO9fjxXkZRZsWZrq+0/eH160dKkpKTx5duBa5Y7c6MfG+kJGn8heML7tuwYUOx5mGBP1166aV6/PHH1axZs5P+2mPGjDluguei3HnnnRo8eLD69OlzwmPh3ofR+mDMoQMAAAAgOAJeyImkPHdmlrWQE+tCCzmIDWPHjlV2drYnBZ2SSklJCVvMKQ0KOgBQHMyxAwAAAPhSixYtTpjQ+WQZM2ZMiZa/9dZby+21mUMHAAAAQHBcd51zQUSrV1+n1avJUSTXvXmdrnuT/CD4GKEDAPkYhQMAgP9t3+51BL536FD55ChWT2W+/Yfw+bHWnnD6aeBkKc38xhR04Cv8nsZJwRsNAAAAIeLj45WTk6O6detS1MFJZ61VTk6O4uPjS/Q8CjoAAAAAgEotKSlJ27dv165du7wOBZVUfHy8kpKSSvQcCjoAKgdG5QAAACCCqlWr6pxzzvE6DKBEKOgAAAAACI6uXb2OwPdq1yZH0XRNIj+IDRR0AAAAAATHX//qdQS+16oVOYrmr33JD2IDBR0AAAAAQLEUPoo9Fs56BQRVFa8DAAAAAIBiu/xy54KIVqy4XCtWkKNILp9+uS6fTn4QfIzQAYCKxq4sAADKT06O1xH43tGjJy9H4c474feuTs5B3kOIDRR0gCKEfkn5/csJAAAAAFA5UNABAAAAgHISxBErAIKJOXQAAAAAAAAChhE6QBjh9qwAAADAB/r08ToC30tMJEfR9DmH/CA2UNCB5yieAAAAoNhGj/Y6At9r3pwcRTP6AvKD2MAhVwAAAAAAAAHDCB0AsYvhXwAAxJ4BA5zruXO9jcPHli1zctS5MzkKZ8CrTn7mXkt+EGwUdAAAAAAEx6FDXkfge8eOkaNoDuWSH8QGCjpACYQO+OD0k4iK0UEAAAAAKhBz6AAAAAAAAAQMI3QAAAAAoAIxcBdARaCgAwAAACA4Lr7Y6wh8r0EDchTNxc3JD2IDBR0AAAAAwXHPPV5H4HtNm5KjaO5JJz+IDcyhAwAAAAAAEDCM0AEQOzhAHQCA2Nezp3OdmellFL62eHFPSVJ6eqancfhVz8k9JUmZGZmexgGUFQUdeILf3QAAAIgF9GsBeIWCDgCcbKE9v9mzvYsDAAAAQGAxhw4AAAAAAEDAUNABAAAAAAAIGA65AgAAABAcV17pdQS+17AhOYrmymTyg9hAQQcAAABAcNx+u9cR+F7jxv7KUbiJo72cRvD2jv7KD1BaHHIFAAAAIDgOHnQuiOjYsYM6dowcRXIw96AO5pIfBB8jdADAS5zxCgCAkhk40LnOzPQ0DD9btszJUXp6preB+NTAV538ZGZkehsIUEaM0AEAAAAAAAgYCjoAAAAAAAABQ0EHAAAAAAAgYJhDB0CwhTttAgAAADxD9ww4OSjoAAAAAAiOjAyvI/C9Ro0yvA7B1zLSMrwOASgXFHQAAAAABAcFnSJR0ImOgg5iBQUdnDQMvQQAAECZ7d7tXCcmehuHjx096uSoWjVyFM7ug05+EquTHwQbBR0AAAAAwTF0qHOdmelpGH62YoWTo/T0TG8D8amh0538ZGZkehsIUEac5QoAAAAAACBgKOgAAAAAAAAEDAUdAAAAAACAgGEOHaAchE74PHu2d3EAAAAAACoHCjoAAAAAguM3v/E6At9r3JgcRfObDuQHsYGCDoDgCR0SBQAAKperrvI6At9r2JAcRXNVCvlBbKCggwrF724AAACUq23bnOtGjbyNw8cOHXJylJBAjsLZts/JT6Na5AfBRkEHAAAAQHBcf71znZnpaRh+tnq1k6P09ExvA/Gp62c6+cnMyPQ2EKCMKOgAgF8wuzYAAACAYqKgAwAAAADFwHQC5SdcLtmfBZRMFa8DAAAAAAAAQMlQ0AEAAAAAAAgYDrkCAAAAEBx33+11BL7XtCk5iuburuQHsYGCDlBKHEONCsUEyQAAhEcnrEgNGpCjaAa1ID+IDRR0AAAAAATHpk3OdYsWFfoyQa4b7d/v5KhmzYrNUVBt2u3kp0Ui+UGwUdABAAAAEBy//rVznZnpaRh+tnatk6P09ExvA/GpX89x8pOZkeltIEAZUdBBuQvy3gwAAAAAAIKAs1wBAAAAAAAEDAUdAAAAAACAgKGgAwAAAAAAEDDMoQMAAAAgOP70J68j8L1mzchRNH/qQX4QGyjoAAAAAAiOvn29jsD36tUjR9H0bUJ+EBs45AoAAABAcGRlORdEtG9flvbty/I6DN/K+jZLWd9meR0GUGaM0AEAAAAQHCNHOteZmV5G4Wvr14+UJKWnZ3oah1+NfG+kJCkzI9PTOICyYoQOAAAAAABAwFDQAQAAAAAACBgKOgAAAAAAAAFDQQcAAAAAACBgmBQZAAAAQHD85S9eR+B7LVuSo2j+0of8IDZQ0AEAAAAQHOnpXkfge3XqkKNo0huRH8QGDrkCAAAAEByLFzsXRPT994v1/ffkKJLF2xZr8Tbyg+BjhA5QzgYN+vn27NnexQEAABCTHnjAuc7M9DQMP9u40clRenqmt4H41AMfOvnJzMj0NhCgjBihAwAAAAAAEDAUdAAAAAAAAAKGQ64AwO84jg8AAABAIRR0AARDaFEDAAAAACo5CjoAAAAAgmP8+ApZbSztO0pOHu91CKVSeBtU1MDk8ReOr5gVAycZBR0AAAAAwZGW5nUEvlerVprXIfha2plpXocAlAsmRQYAAAAQHPPmORdEtGvXPO3aRY4imfflPM37kvwg+BihAwAAACA4HnnEue7b19s4fOzzz50c1atHjsJ5ZIGTn75NyA+CjYIOAP+KpYPZAQAAAKAcUdBBueB3NwAAAAAAJw9z6AAAAAAAAAQMI3QAAAAAAL4T7iiAijqVORBEFHSAChT6JcSXDwAAQDn4xz+8jsD3UlPJUTT/uJj8IDZQ0AGAICm8q4pKIQCgsmnRwusIfK9mTXIUTYtE8oPYwBw6AAAAAIJj9mx2aBRh587Z2rmTHEUye9Nszd5EfhB8jNABAAAAEBxPPOFcc5rViL74wslRgwbkKJwnljj5GdSC/CDYGKEDAAAAAAAQMBR0AAAAAAAAAoZDrgAAAABUKhytBSAWMEIHAAAAAAAgYBihA8Bf2GUGAACieeUVryPwvbZtyVE0r1xKfhAbKOgAAAAACI5GjbyOwPcSEshRNI1qkR/EBgo6wEkSOvBk9mzv4gAAAAi0adOc66uu8jYOH9uxw8lRw4bkKJxp65z8XJVCfhBsFHQAIMioFAIAKpu//925pqAT0datTo4o6IT39xVOfijoIOiYFBkAAAAAACBgKOgAAAAAAAAEDAUdAAAAAACAgKGgAwAAAAAAEDBMigwAAAAgON54w+sIfK9DB3IUzRtXkh/EBgo6AAAAAIIjMdHrCHyvWjVyFE1idfKD2EBBByUSeoZkAAAA4KSbPNm5zsgo9lMqWx9227bJkqRGjTI8jcOvJmdNliRlpGV4GgdQVsyhAwAAACA4Jk/+uaiDsLZtm1xQ1MGJJmdNLijqAEHGCB0AAAAAQCCEG201e/bJjwPwA0boAAAAAAAABAwFHQAAAAAAgIChoAMAAAAAABAwzKEDeCD02F+O+QUAACiBd9/1OgLf69y5cuWo8Lw6RfWv3722cuUHsYuCDgDECiqFAIDKoHp1ryPwvbg4chRN9arkB7GBQ64AAAAABMeECc4FEW3dOkFbt5KjSCZ8MkETPiE/CD5G6AAAAAAIjunTnevbbw/7cLjTWlc2O3Y4OWrcOHyOKrvp65383N6R/CDYGKEDAAAAAAAQMIzQAeA9dqUBAAAAQIkwQgcAAAAAACBgKOgAAAAAAAAEDIdcAQAAAAiOzEyvI/C99PRMr0PwtcyMTK9DAMoFI3QAAAAAAAAChhE6AAAAAIJj3Djn+p57vI3Dx774wslR06aVM0fhzrcxe/bPt8ctdvJzT3rlzA9iByN0AAAAAATHnDnOBRHt3DlHO3eSo0jmbJ6jOZvJD4KPgg4AAAAAAEDAUNABAAAAAAAIGAo6AAAAAAAAAcOkyChSuEnFAAAAAE8kJHgdge/FxZGjaBKqkh/EBgo6AAAAAIJj7lyvI/C9zp3JUTRzryU/iA0UdACPhY6ACj2dYkxj2BcAACgHdCkAVGYUdAAgFlXKSiEAoFJ4+GHnevRob+Pwsc2bnRw1b06Ownn4P05+Rl9AfhBsTIoMAAAAIDg+/NC5IKLduz/U7t3kKJIPv/pQH35FfhB8FHQAAAAAAAAChoIOAAAAAABAwFDQAQAAAAAACBgmRQYAAAAQHHXreh2B71WrRo6iqVud/CA2UNABAAAAEBwzZngdge916ECOoplxJflBbKCgA+DkCD2NNgAAAACgTCjoAECsCy2mzZ7tXRwAAJSH++93rv/6V2/j8LENG5wctWpFjsK5f56Tn7/2JT8INgo6AAAAAIJjyRKvI/C9PXvIUTRLtpMfxAYKOgAAAACAmBY6YPnTFlLrFO9iAcoLpy0HAAAAAAAIGAo6AAAAAAAAAcMhVwiLExIBAADAl5KSvI7A9xISyFE0CUeTlHS611EAZUdBBwAAAEBw/POfXkfge23bkqNo2n71T/3zGa+jAMqOQ64AAAAAAAAChoIOAAAAgOAYOdK5IKL160dq/fqRXofhW+sbjdTI90Z6HQZQZhxyBfhI6NxFs2d7FwcAAIBvZWV5HYHv7duX5XUIvravepayvvU6CqDsGKEDAAAAAAAQMBR0AAAAAAAAAoZDrgBUnNBjyAAAAAAA5YaCDgAAAIBAGDRI+u1XzSVJz7HfKKKaNZt7HYKv1TzcXM3reh0FUHYUdAAAAAAExnOpz3sdgu+lkqOoUr9+Xs9TEEQMYA4dAAAAAACAgKGgAwAAACAwfrv2Nv127W1eh+Fra9feprXkKKK1v7xNt80mPwg+DrkCfCp0PuHZs72LAwAAwE/O3r/Z6xB8bz85imp//GZtzvE6CqDsKOgAQGVCpRAAAACICRxyBQAAAAAAEDAUdAAAAAAAAAKGQ65QYBCn7gMAAIDPfVkrzesQfK8WOYqq1sE0fbNcGjTj+Ps5Gh1BQ0EHAAAAQGC8mDze6xB8L5kcRZW8bbzXIQDlgkOuAAAAAAAAAoYROkAAcGIiAAAAx12rr5MkPdn2nx5H4l+r3Ry1JUdhrT7Hzc9X5AfBRkEHAAAAQGAkHtrudQi+d4gcRXWoGvlBbKCgAwAAAACo9MKdJIbR8fAzCjqVHGe2AgAAAAAgeJgUGQAAAAAAIGAYoQMAAADAdyKNJN9Yu+vJDSSAapOjqGrvJz+IDRR0AAAAAATGy63+6nUIvteKHEXV6r/kB7GBQ64AAAAAAAAChoIOAAAAgMC4f8Xlun/F5V6H4WsrVlyuFeQoohVNL9eKpuQHwcchVwAAAAAC47SjOV6H4HtHyVFUR08hP4gNFHQAAAAAAAij8OTcs2d7EwcQDgUdAOUr0ikp4D+h24reCQDAY3QhAKBkmEMHAAAAAAAgYBihAwAAACAw1iT28ToE30skR1El/kB+EBso6AAAAAAIjGnNR3sdgu81J0dRNc8mP4gNFHQAlB0HvQcfM/4BAAAAgcIcOgAAAAACY8yyARqzbIDXYfjasmUDtIwcRbSs2QAta0Z+EHyM0KmEGEwBAACAoKp27JDXIfjeMXIU1bEq5AexgRE6AAAAAAAAAUNBBwAAAAAAIGAo6AAAAAAAAAQMc+gAAAAACIxPGlzsdQi+14AcRdVgL/lBbKCgAwAAACAwZja9x+sQfK8pOYqq6U7yg9hAQQcImNCzlM2e7V0cAAAAAADvMIcOAAAAgMD4y+Ke+svinl6H4WuLF/fUYnIU0eIWPbW4RU+vwwDKjIIOAAAAAABAwHDIFRBgHH6FCsObCwAAAPA1RugAAAAAAAAEDAUdAAAAAACAgOGQKwAAAACBsbDhlV6H4HsNyVFUDb8nP4gNFHQAAAAABMa7jW/3OgTfa0yOomq8i/wgNlDQAVA8oZPkSkyUCwAASq1wt6IkTj12UJJ0JK56OUUTe465OYojR2Edq+Lm56eS5yfce5duMbxCQQdA6ZSlJwYAAFBKDy0bKEl6ID3T20B8bJmbo3RyFNayZm5+NmV6GwhQRhR0Kgl+ewMAAAAAEDso6AAAogutCDOmGAAAAPAFCjpAjOA3NwAAAABUHlW8DgAAAAAAAAAlwwgdAEDxMRQMAFBC5T2X44eNMsp3hTGoETmKqtHuDK9DAMoFBR0AAAAAgUFBp2gUdKJrlJPhdQhAueCQKwAAAACBcfrR3Tr96G6vw/C1o0d36yg5iujoKbt19BTyg+BjhA4AAACAwLhvxVBJ0gPpmd4G4mMr3Bylk6OwVjR187Mp09tAgDJihA4AAAAAAEDAUNABAAAAAAAIGA65imHlfUYBAAAAAMDxCv/u4kSgOFkYoQMAAAAAABAwjNABAAAAEBjvNv6N1yH4XmNyFFXj78gPYgMFHQAAAACBsbDhVV6H4HsNyVFUDfeQH8QGDrkCAAAAEBiJh7Yp8dA2r8PwtUOHtukQOYroUNVtOlSV/CD4GKEDAAAAIDDuWn29JOmB9ExvA/Gx1W6O0slRWKubuPnZlOltIEAZUdABAAAAUG440yoAnBwccgUAAAAAABAwjNCJMewRgXT8+2D2bO/iAAAAAABUDAo6ACKjQggAAAAAvkRBBwAAAECpeLHvZ2bTu0/+iwZMU3IUVdNvyQ9iAwUdAAAAAIHxSQNGEBelATmKqsG+is1PuEIn0yCgIjApMgAAAIDAOHv/Jp29f5PXYfja/v2btJ8cRbT/1E3afyr5QfAxQgcAAABAYPx27a8lSQ+kZ3obiI+tdXOUTo7CWtvYzc+mTG8DAcqIgg4AoHQ4nRoAAADgGQo6QIzjNzcAAADgLebVQUWgoBMDOLM0AM9ROQSASoF+JwD4B5MiAwAAAAAABAwjdAAcj11vAADAx6Y1+5PXIfheM3IUVbMd5AexgYIOAAAAgMBYU6+v1yH4Xj1yFFW9H8kPYgOHXAEAAAAIjHP2ZemcfVleh+Fr+/ZlaR85imhfQpb2JWR5HQZQZozQAQAAABAYt64fKUl6ID3T0zj8bL2bo3RyFNb6X4yUJKVvyvQ0jsIzHXBeCZQUBR2gEol4IiLmzQEAAACAQOGQKwAAAAAAgIChoAMAAAAAABAwFHQAAAAAAAAChjl0AoopTwAAAFAZvdzyL16H4HstyVFULbeTH8QGCjoAAAAAAmNjnXSvQ/C9OuQoqjoHyA9iAwUdoLJimBcAAIjCr12Flt8vlkRhJ5rv3RxR2Anv+xpufijsIOAo6AAAAAAIjBs2PiBJeiA909tAfGyjm6N0chTWxiQ3P5syvQ0EKCMKOgAAAAAAeCzcqLjZs09+HAgOCjpAJbV8+c+3O3XyLg4AAOA9vx5eBQCIjIIOAKB8Ff5VwK4lAAAAoNxV8ToAAAAAAAAAlAwjdIBKZPRyxlPDA6EjdhitAwAooxeSx3sdgu8lk6Ookr8Z73UIQLmgoAMAAAAgML6qleZ1CL5XixxFVetQmtchAOWCQ64AAAAABEabXfPUZtc8r8PwtV275mkXOYpo12nztOs08oPgY4ROgHD2AVQUzngFAACC4qrPH5EkranX1+NI/OtzN0f1yFFYnzd087PJ//nhXBOIhhE6AAAAAAAAAUNBBwAAAAAAIGA45AoAAACoZDiUHwCCjxE6AAAAAAAAAcMIHQAAAACB8VzqP7wOwfdSyVFUqVuDm59wo+uYKLnyoqDjcwyHxcnGGa8AAICf/bdmC69D8L2a5CiqmkfID2IDBR0gho1eTkUQAADElo47neEInzSgnxPJTjdHDchRWDtrufnZR34QbBR0fIhROQBiVugHHOODAeCkiLW+5aVfPCGJgk40X7g5oqAT3hdnuvmhoIOAo6ADIKLQw68kDsFCOaO4AwAAAJQaZ7kCAAAAAAAIGEboADGGeXMAAAAAIPZR0PGJWDu2GbGJM2ABAAAA/sKpzCsvCjoAAO8xnw4AoJiebPuK1yH4XltyFFXbL8kPYgMFHSAGcJgVAACQKseo790JjbwOwfcSyFFUCbnkB7GBgg4AAAAQUJWhgFNYtx3TJEkLG17lcST+tcPNUUNyFNaO2m5+9pAfBBsFHQAAAACBMXDr3yVR0Ilmq5sjCjrhba3v5ieGCzqFi70c0R6bKOgAKBUmSEaFYT4dAAAAoEgUdDxUGYfIAgAAAACAsqOgAwSUnyZCZrQOKgzjhQEAAICwKOgAAAAAAcDobgClFe7zg/1kwUdB5yTjixgAAAAovbEd3vA6BN/rQI6i6vAF+ZEo8sQCCjpAgPjpMCsAAAAv/FAt0esQfK8aOYqqWh75QWygoHMSMCoHAMoJZ8ACUInQhwyvz7bJkqQPG2V4GoefbXNz1IgchbWt7mRJUqOcDE/jAMqKgg7gc4zKAQAg9lG8KT4KOkWjoBPdtsTJkijohMP5KIKFgg6AcsUZr3DSMFoHAAAAlRgFHcAHGIUDAAAAACgJCjrliKGyAAAAAADgZKCgA6DCcPgVAAAAEFzFGbTAke/eoaAD4KQILe6EotCDcsF8OgBQafxv53e9DsH3OpOjqDp/Tn4QGyjolBCHVQEAAADeORJX3esQfC+OHEUV9xP5QWygoAOcREx+fCIOy0K5izRah1E8AE4Cdv5VvIFbJ0iS3m18u8eR+NdWN0eNyVFYW+u5+dlFfhBsFHQA+AbFHZQ7flkBKEd8pPhDtx3TJVHQiWaHmyMKOuHtqOPmh4JOuQj32ci+s5ODgk6I4uzUBXByUNxBhYr2wU4PBIDo/wFAWVDkOTko6ETAlzjKgkOrgADj0CwAAIByV16/seme/azSFHQYfYPSCi3OPNwp8qcHRZyKE+kMWZEwogcVItIXBr0KwPcK//vybwsAiAWBKuiUdJ7LSH1vijiVS7RCS2iBpjgFGYo2weBVAYjDxGJQcb4wivulUlETNDOiCDHgZPfN6AsCQHCd7M/wwt0rPx1OZqy1ZV+JMbskfR3moURJu8v8AsFCm2NfZWuvRJsrC9pcOZS2zb+01tYr72BQNlH6YH5VGf/nQlXm9lfmtkuVu/2Vue1S5W5/ZW67VH7tj9gHK5eCTiTGmBXW2g4V9gI+RJtjX2Vrr0SbKwvaXDlUxjbDPyr7+68yt78yt12q3O2vzG2XKnf7K3PbpZPT/ioVuXIAAAAAAACUPwo6AAAAAAAAAVPRBZ3nK3j9fkSbY19la69EmysL2lw5VMY2wz8q+/uvMre/Mrddqtztr8xtlyp3+ytz26WT0P4KnUMHAAAAAAAA5Y9DrgAAAAAAAAKmWAUdY8ydxph1xpj1xpiRYR7vaYzZZ4zJci8Phjx2oTFmkzFmizHmvpD7XzXGrDXG/CXkvj8ZYy4pW5NKxxjzkjHmO2PMupD76hhjPjDGfO5e147w3BvdZT43xtwYcn97Y8ynbtufMcYY9/7H3La/HLLsdeFyW5HK2OZjIdv77ZD7zzHGLHPbPM0YU829/w73PfRuyH3djDFPVXQ7C8Udrs1XuO/tn4wxEWchj/Je9m2by9jere77N8sYsyLk/rDvEWPM5e56PzbG1HXva2qMmVaRbQwTd7g2/z9jzEb3/26mMeaMCM8N3DZ2X7MsbY6l7fyw294sY8z7xpiGEZ4bS5/ZxW1zID+z4X+RPivCLBcz78FI3xUhj5/qtmWL27bGIY/d796/yRjT372vnjFmodveS0KWnRXpf9pLxWh/hjFmV8j2viXksRM+f918vee2//aQZZ83xrQ7Oa0qnnCfw4UeN+73xxb3s7ldyGOBbrtUrPYH/jdhJMaYRsaY+caYz4zTD7ozzDIxuf2L2fZY3vbxxpjlxpg1bvv/N8wy3n3uW2ujXiSlSFonqbqkUyTNk3RuoWV6SpoT5rlxkr6Q1ERSNUlrJJ0nKVXSi+4yH0iqJeksSbOLiqeiLpJ6SGonaV3IfY9Lus+9fZ+kx8I8r46kL93r2u7t2u5jyyV1kWQkzZU0wG3rB+7jL0pqLSlB0oeSqgahze5j+yPcP13S1e7tiZJ+495eKqeA+CdJg9yc/FtSHR+0uZWkFpIyJXWI8Lyw72W/t7m07XWX2yopMcz9Yd8j7vqqS7pO0h3ufVMlNfPBNv6VpFPc249F+F8O5DYuS5tjcDufHnL795ImhnlerH1mF9lm97FAfmZz8f8l0mdFmOVi4j0Y7bsiZJnb8/8XJV0taZp7+zx3+VMlneOuJ879373O/WzNdJcdJGmM19u3lO3PkPS3MM8N+/krabC7natIWuIu20bSJK/bG6YNJ3wOF3p8oPv9Ydzvk2Wx0vZitr+nAv6bMErbz5LUzr19mqTNYd77Mbn9i9n2WN72RlJN93ZVScskdSm0jGef+8UZodPKfTMetNbmSfqPpMuK8TxJ6iRpi7X2S2vtUUmvSRoiKVdSgjGmipuUY5L+LOmhYq633FlrF0j6vtDdQyRNcW9PkXRJmKf2l9PZ/95au0fOG/JCY8xZcjraS62zhV52n/+TpKrGGCNnA+ZKukfSs9ba3PJtVXRlaHNYbpt6S3ojzPONnG2d3+brJM211hZ+/QoVrs3W2g3W2k1FPDXse9nvbS5De6OJ9B75Sc6HVXVJucaY7pK+tdZ+XobXKrEIbX7f/fySnB8JSWGeGshtLJWpzdEEcTv/EPJnDUnhJomLqc/sYrY5rCC8txEIMd1vCCNS3zZUaE7ekNTHbesQSa9Za49Ya7+StMVdX66cdp4q6Zgx5hRJI+UUy/ymOO2PJOznr35uf1U5212SHpY0ulwjLwcR+s6hhkh62TqWSjrD/X4JfNulYrU/ksD8JozEWpttrV3l3v5R0gZJZxdaLCa3fzHbHkksbHtrrd3v/lnVvRTub3n2uV+cgs46Sd2NMXWNMdXlVB4bhVmuqzsMaa4xJtm972xJ20KW2S7pbGvtBkm7JK2SNFvSuZKq5L9RfKSBtTbbvf2tpAZhlgnbRveyvfD97j/Bu5JWS8qWtE9SZ2vtW+UbeqkVp82SFG+MWWGMWRoyTKyupL0hPyLzcyFJf5Pzo/IXkhZJuknSc+UdfAWKtJ1juc1W0vvGmJXGmNtC7o/0HvmrnBF8g+SM2Bgt50vJb4bL2XtSWCxv40htlmJsOxtjHjXGbJN0raQHwywSa5/ZxWmzVPk+s3HyVLZ+Q6TPkLDLuG3bJ6etkZ77Lzmd/g8k/UXOnt5XrLUHKyD+sipO+yXpcvdQijeMMfm/GyI99wNJjeVs72eMMYMlrbLW7ijv4E+CaN8xsd72fLH4m/A47uE0beWM1AgV89s/StulGN72xpg4Y0yWpO/kFOcibvuT/bl/SlELWGs3GGMek/S+pAOSsuRU0EKtkvRLa+1+Y8xASW9JalbEekfm3zbGzJb0a2PMKDnDzD6w1r5Q/GZUPGutNcYUe89nEet6XG71zRjzoqQHjXN88a8krbXWPlIer1NWRbT5l9ba/xpjmkj6yBjzqZw3bqR1vSLpFUlyj6l8RtIAY8wNct7kd1trfyrfFngrBtrczd3G9SV9YIzZ6O6ZKRD6HrHWfiDnQ0luG9+V1NwYc4+kPZLu9Lpz6n7G5El6tTzWF4RtXIw2x9R2ttaOkjTKGHO/pN+pHPb0+P0zu5ht5jMbpWaMmSfpzDAPjQr9g35D6Vhr90m6SJKMMwfRfZIuNca8IOewjCestUs8DLGkZkuaaq09Yoz5tZy91r0jLez++LlGkowxVeUcWjfEGPOknILey9batyM9P8hisO0x/5vQGFNT0gxJI+3xo2RLLGjbv4i2x/S2t9Yek5RmnDkpZxpjUqy1YeeSKub6yu1zv1iTIltrJ1lr21tre8jpsG8u9PgP+cOQrLXvyhmenijpvzp+NE+Se18BY8wQSSsl1ZTU1Fp7paSh7mggr+10h8nJvf4uzDKR2vhfHX+IQ7i2t5UzvG6TpCvctjc1xkR981ew4rRZ1tr/utdfyplTo62kHDlDC/MLheHa3FBSJ3fv9t2SrpK0V1Kfcm5HeYu0nWO2zSHb+DtJM+UMD5SKeI+4/7sZcvak/q+kGyUtlDN6wDPGmAxJF0u61lob7gdHzG3jYrQ55rZziFclXR7m/lj7zA4Vqc2V8TMb5cha29damxLmMkuVr99QZN82dBm3bbXktLU4zx0t6VFJw+R8pt4oaUz5hF4uimyDtTbHWnvE/fNFSe2L+1w5e6lfljP/yD452/vucon85Ij2HRPrbY/l34SSCoouMyS9aq19M8wiMbv9i2p7rG/7fNbavZLmyzlkLpRnn/vFPctVfff6F3Lmz/lXocfPNKbgbCCd3PXmSPpEUjPjnMGgmpwJgkLPbFBVPx8rlqCfj0WLkzNpktfelpNQudezwizzb0m/MsbUdqtrv5L0b+sMP/7BGNPFzc0NYZ6ff4xkVTltlpz5Grx88xbZZretp7q3EyWdL+kz9wfjfElDozz/Yf18SED+Nve6zcUR9r0cq202xtQwxpyWf1vO+zq/Cl3Ue+QPkp6xzvwivmivMeZCSX+UNDjK6JGY2sbFaXMMbufQwsoQSRvDLBZTn9nFaXMl/czGyVPZ+g1R+7au0JwMlfSR29a3JV1tnLOhnCNn7/Xy/Ce5/89J1tpMOe37SU57EyqwPSVVZPvzC3yuwXLm25AifP6GPK+2nJ0QL8u/7S/K25JuMI4ukva53y+Voe2x/Jswf86vSZI2WGufjLBYTG7/4rQ9xrd9PeOeLdYYkyCpn07sb3n3uW+LN7Pzx5I+kzNDcx/3vhGSRri3fydpvfv4UknpIc8dKGdEzxeSRhVa70hJGfbn2aOnSvpUEc6QUJEX97Wz5UxQtF3SzXKOe/tQ0udy5oyo4y7bQe6s3O7fw+VMcLRF0k0h93eQ8+PoCznHgpuQxy5RyCzWksa5bX/V722WlO7Gusa9vjlknU3kvEm3SHpd0qkhj7VVyKzt7vZfL+m90OU8aPOl7u0jknbK+XEnSQ0lvVvUe9nPbS5te902rXEv6wu1N+x7JGQd74T8fYX7/EWS6nm4jbfIGaKf5V4mhsQb6G1cljbH4HaeIeczd62cIf9nu8vG8md2kW1WgD+zufj/EumzIpbfgwrzXSFnMs/B7u14ty1b3LY1CXnuKPd5myQNKLTe6XLPGCipvqTFbnsv93o7l7D9f9XPvwvmS2oZ8tywn7/uY09J6hmSw/fd9dzhdZtDYgz3ORz6m8jIGbn6hfte7xDy3EC3vZjtD/xvwiht7ybnh/Za/dy3GlgZtn8x2x7L2z5VzlyKa+X0uR507/fF575xnwwAAAAAAICAKNYhVwAAAAAAAPAPCjoAAAAAAAABQ0EHAAAAAAAgYCjoAAAAAAAABAwFHQAAAAAAgIChoAMAAAAAABAwFHQAAAAAAAAChoIOAAAAAABAwPx/eytvQbrPndwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 2, figsize=(20, 6));\n",
    "axs[0].hist(control_samples, density=True, bins=100, alpha=0.7, color='b', label='Control');\n",
    "axs[0].hist(test_samples, density=True, bins=100, alpha=0.7, color='r', label='Treatment');\n",
    "axs[0].xaxis.set_major_formatter(ticker.PercentFormatter(xmax=1));\n",
    "axs[0].set_title('Session Conversion Rate');\n",
    "axs[0].get_yaxis().set_visible(False);\n",
    "axs[0].legend();\n",
    "\n",
    "\n",
    "diffs = test_samples - control_samples\n",
    "p10_diff = np.percentile(diffs, 10)\n",
    "p90_diff = np.percentile(diffs, 90)\n",
    "mean_diff = np.mean(diffs)\n",
    "axs[1].hist(diffs, density=True, bins=100, alpha=0.7, color='b');\n",
    "axs[1].axvline(p10_diff, color='r', linestyle='dashed', label=f'10th Percentile ({p10_diff*100*100:0.0f} bps)');\n",
    "axs[1].axvline(mean_diff, color='k', linestyle='dashed', label=f'Mean ({mean_diff*100*100:0.0f} bps)');\n",
    "axs[1].axvline(p90_diff, color='g', linestyle='dashed', label=f'90th Percentile ({p90_diff*100*100:0.0f} bps)');\n",
    "axs[1].xaxis.set_major_formatter(ticker.PercentFormatter(xmax=1));\n",
    "axs[1].set_title('Difference in Session Conversion Rate');\n",
    "axs[1].get_yaxis().set_visible(False);\n",
    "axs[1].legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adjusting for nonindependence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T21:09:10.512054Z",
     "iopub.status.busy": "2021-01-31T21:09:10.511642Z",
     "iopub.status.idle": "2021-01-31T21:09:11.167710Z",
     "shell.execute_reply": "2021-01-31T21:09:11.167085Z",
     "shell.execute_reply.started": "2021-01-31T21:09:10.512031Z"
    }
   },
   "outputs": [],
   "source": [
    "data = {\n",
    "    'user': [],\n",
    "    'session_id': [],\n",
    "    'assignment': [],\n",
    "    'session_converted': []\n",
    "}\n",
    "\n",
    "# create a copy of the session conversion rates for each user\n",
    "# we'll use this to keep track of each user's session conversion rate\n",
    "# throughout the simulation in case it changes\n",
    "conversion_rates = baseline_conversion_rates.copy()\n",
    "\n",
    "# Simulate all sessions for each user\n",
    "for user_id, num_sessions in enumerate(sessions_per_user):\n",
    "    for session_id in range(1, num_sessions+1):\n",
    "        # randomly assign session to control (0) or test (1)\n",
    "        assignment = np.random.randint(0, 2)\n",
    "        \n",
    "        if assignment == 1:\n",
    "            # increase user's conversion rate permanently\n",
    "            conversion_rates[user_id] = baseline_conversion_rates[user_id] + conversion_uplifts[user_id]\n",
    "            new_conversion_rate = conversion_rates[user_id]\n",
    "        else:\n",
    "            new_conversion_rate = conversion_rates[user_id]\n",
    "        \n",
    "        # see if the session converted\n",
    "        session_converted = np.random.choice([0, 1], p=[1-new_conversion_rate, new_conversion_rate])\n",
    "        \n",
    "        # record the results\n",
    "        data['user'].append(user_id)\n",
    "        data['session_id'].append(f\"{user_id}-{session_id}\")\n",
    "        data['assignment'].append(assignment)\n",
    "        data['session_converted'].append(session_converted)\n",
    "        \n",
    "df = pd.DataFrame(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Local Dask Cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T21:09:11.170059Z",
     "iopub.status.busy": "2021-01-31T21:09:11.169879Z",
     "iopub.status.idle": "2021-01-31T21:09:12.394881Z",
     "shell.execute_reply": "2021-01-31T21:09:12.394096Z",
     "shell.execute_reply.started": "2021-01-31T21:09:11.170037Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table style=\"border: 2px solid white;\">\n",
       "<tr>\n",
       "<td style=\"vertical-align: top; border: 0px solid white\">\n",
       "<h3 style=\"text-align: left;\">Client</h3>\n",
       "<ul style=\"text-align: left; list-style: none; margin: 0; padding: 0;\">\n",
       "  <li><b>Scheduler: </b>tcp://127.0.0.1:59867</li>\n",
       "  <li><b>Dashboard: </b><a href='http://127.0.0.1:8787/status' target='_blank'>http://127.0.0.1:8787/status</a></li>\n",
       "</ul>\n",
       "</td>\n",
       "<td style=\"vertical-align: top; border: 0px solid white\">\n",
       "<h3 style=\"text-align: left;\">Cluster</h3>\n",
       "<ul style=\"text-align: left; list-style:none; margin: 0; padding: 0;\">\n",
       "  <li><b>Workers: </b>12</li>\n",
       "  <li><b>Cores: </b>12</li>\n",
       "  <li><b>Memory: </b>34.36 GB</li>\n",
       "</ul>\n",
       "</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Client: 'tcp://127.0.0.1:59867' processes=12 threads=12, memory=34.36 GB>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cluster = LocalCluster(\n",
    "    n_workers=12, \n",
    "    threads_per_worker=1, \n",
    "    processes=True, \n",
    "    silence_logs=logging.ERROR\n",
    ")\n",
    "client = Client(cluster) # tcp://127.0.0.1:59867\n",
    "client\n",
    "\n",
    "# client.shutdown()\n",
    "# client.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run Simulation Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T23:43:54.717762Z",
     "iopub.status.busy": "2021-01-31T23:43:54.709514Z",
     "iopub.status.idle": "2021-01-31T23:43:54.911382Z",
     "shell.execute_reply": "2021-01-31T23:43:54.892502Z",
     "shell.execute_reply.started": "2021-01-31T23:43:54.717716Z"
    }
   },
   "outputs": [],
   "source": [
    "def calculate_metrics(data, metric):\n",
    "\n",
    "    df = pd.DataFrame(data)\n",
    "    num_samples = 50000\n",
    "    \n",
    "    if metric == 'user':\n",
    "        df = df.groupby(['user', 'assignment']).agg(converted=('session_converted', 'max')).reset_index()\n",
    "    elif metric == 'session':\n",
    "        df['converted'] = df['session_converted']\n",
    "    else:\n",
    "        raise NotImplementedError(\"Metric must be 'user' or 'session'\")    \n",
    "\n",
    "    control_converted = df[df.assignment==0].converted.sum()\n",
    "    control_total = df[df.assignment==0].converted.count()\n",
    "    control_samples = np.random.beta(control_converted, control_total - control_converted, size=num_samples)\n",
    "\n",
    "    test_converted = df[df.assignment==1].converted.sum()\n",
    "    test_total = df[df.assignment==1].converted.count()\n",
    "    test_samples = np.random.beta(test_converted, test_total - test_converted, size=num_samples)\n",
    "\n",
    "    # how often does test convert higher than the control?\n",
    "    test_gt_control = (test_samples > control_samples).mean()\n",
    "    # in what range do we think the true conversion rates lie?\n",
    "    test_interval = [np.percentile(test_samples, 10), np.percentile(test_samples, 90)]\n",
    "    control_interval = [np.percentile(control_samples, 10), np.percentile(control_samples, 90)]\n",
    "    \n",
    "    # how much better is test than control?\n",
    "    diffs = test_samples - control_samples\n",
    "    diff_interval = [np.percentile(diffs, 10), np.percentile(diffs, 90)]\n",
    "    mean_diff = np.mean(diffs)\n",
    "    \n",
    "    return {\n",
    "        'test_gt_control': test_gt_control,\n",
    "        'test_interval': test_interval,\n",
    "        'control_interval': control_interval,\n",
    "        'diff_interval': diff_interval,\n",
    "        'mean_diff': mean_diff\n",
    "        \n",
    "    }\n",
    "        \n",
    "def run_simulation(\n",
    "        num_users=30000, \n",
    "        baseline_conversion=0.1, \n",
    "        relative_uplift=0.1, \n",
    "        is_aa_test=False,\n",
    "        randomization_unit='session',\n",
    "        persist_treatment_effect=False,\n",
    "        geometric_p=0.5,\n",
    "        conversion_dependent_on_spu=False,\n",
    "        conversion_spu_func=None\n",
    "    ):\n",
    "\n",
    "    if randomization_unit not in ['session', 'user']:\n",
    "        raise NotImplementedError(\"randomization_unit must be 'session' or 'user'\")\n",
    "    \n",
    "    # sessions per user\n",
    "    sessions_per_user = np.random.geometric(geometric_p, size=num_users)\n",
    "    \n",
    "    # each user has some baseline conversion rate\n",
    "    if conversion_dependent_on_spu:\n",
    "        _baseline_conversion_rates = conversion_spu_func(baseline_conversion, sessions_per_user)\n",
    "        a = _baseline_conversion_rates*1000\n",
    "        b = 1000 - _baseline_conversion_rates*1000\n",
    "        baseline_conversion_rates = np.random.beta(a, b, size=num_users)\n",
    "    else:\n",
    "        a = baseline_conversion*1000\n",
    "        b = 1000 - a\n",
    "        baseline_conversion_rates = np.random.beta(a, b, size=num_users)\n",
    "\n",
    "    if is_aa_test:\n",
    "        conversion_uplifts = np.zeros(shape=num_users)\n",
    "    else:\n",
    "        # each user will experience a conversion uplift of 10% relative, so 0.2% or 20bps (absolute)\n",
    "        # larger sample size to emphasize confidence in 20 bps effect size for the sake of demonstration\n",
    "        a = baseline_conversion*relative_uplift*10000\n",
    "        b = 10000 - a\n",
    "        conversion_uplifts = np.random.beta(a, b, size=num_users)\n",
    "\n",
    "    # session level randomization\n",
    "    data = {\n",
    "        'user': [],\n",
    "        'session_id': [],\n",
    "        'assignment': [],\n",
    "        'session_converted': []\n",
    "    }\n",
    "\n",
    "    conversion_rates = baseline_conversion_rates.copy()\n",
    "\n",
    "    for user_id, num_sessions in enumerate(sessions_per_user):\n",
    "        if randomization_unit == 'user':\n",
    "            assignment = np.random.randint(0, 2)\n",
    "            \n",
    "        for session_id in range(1, num_sessions+1):\n",
    "            if randomization_unit == 'session':\n",
    "                assignment = np.random.randint(0, 2)\n",
    "            \n",
    "            if assignment == 1 and persist_treatment_effect:\n",
    "                # increase user conversion rate permanently\n",
    "                conversion_rates[user_id] = baseline_conversion_rates[user_id] + conversion_uplifts[user_id]\n",
    "                user_conversion_rate = conversion_rates[user_id]\n",
    "            elif assignment == 1:\n",
    "                # temporarily increase conversion rate for given session\n",
    "                user_conversion_rate = baseline_conversion_rates[user_id] + conversion_uplifts[user_id]\n",
    "            else:\n",
    "                user_conversion_rate = conversion_rates[user_id]\n",
    "            \n",
    "            session_converted = np.random.choice([0, 1], p=[1-user_conversion_rate, user_conversion_rate])\n",
    "            data['user'].append(user_id)\n",
    "            data['session_id'].append(f\"{user_id}-{session_id}\")\n",
    "            data['assignment'].append(assignment)\n",
    "            data['session_converted'].append(session_converted)\n",
    "\n",
    "    # now do bayesian analysis with our observered results\n",
    "    return {\n",
    "        'session_conversion': calculate_metrics(data, metric='session'),\n",
    "        'user_conversion': calculate_metrics(data, metric='user'),\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T23:44:45.907061Z",
     "iopub.status.busy": "2021-01-31T23:44:45.900171Z",
     "iopub.status.idle": "2021-01-31T23:44:51.987210Z",
     "shell.execute_reply": "2021-01-31T23:44:51.969493Z",
     "shell.execute_reply.started": "2021-01-31T23:44:45.907020Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'session_conversion': {'test_gt_control': 1.0,\n",
       "  'test_interval': [0.0902155555788382, 0.09449700302468118],\n",
       "  'control_interval': [0.07602173045769126, 0.08000847063171752],\n",
       "  'diff_interval': [0.011415089259289742, 0.017256651469247014],\n",
       "  'mean_diff': 0.014338944066677486},\n",
       " 'user_conversion': {'test_gt_control': 1.0,\n",
       "  'test_interval': [0.16726041753239856, 0.17516143925513078],\n",
       "  'control_interval': [0.14204752243190533, 0.1494356955346826],\n",
       "  'diff_interval': [0.020084265895938532, 0.030893171551643225],\n",
       "  'mean_diff': 0.025462020520587183}}"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run_simulation(\n",
    "    num_users=30000, \n",
    "    baseline_conversion=0.1, \n",
    "    relative_uplift=0.1, \n",
    "    is_aa_test=False,\n",
    "    randomization_unit='user',\n",
    "    persist_treatment_effect=False,\n",
    "    geometric_p=0.5,\n",
    "    conversion_dependent_on_spu=True,\n",
    "    conversion_spu_func=lambda x,y: x*2/3 + x*1/3*np.exp(-0.4*y)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-24T21:35:54.095657Z",
     "iopub.status.busy": "2021-01-24T21:35:54.095197Z",
     "iopub.status.idle": "2021-01-24T21:38:36.389203Z",
     "shell.execute_reply": "2021-01-24T21:38:36.386747Z",
     "shell.execute_reply.started": "2021-01-24T21:35:54.095622Z"
    }
   },
   "source": [
    "## A/B Test - Session Level Randomization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T19:19:18.942604Z",
     "iopub.status.busy": "2021-01-31T19:19:18.929945Z",
     "iopub.status.idle": "2021-01-31T19:20:35.868713Z",
     "shell.execute_reply": "2021-01-31T19:20:35.861340Z",
     "shell.execute_reply.started": "2021-01-31T19:19:18.942562Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Detected positive effect 100.0% of the time\n",
      "Average effect size detected: 201 bps\n",
      "\n"
     ]
    }
   ],
   "source": [
    "delayed_results = []\n",
    "\n",
    "NUM_EXPERIMENTS = 250\n",
    "for iteration in range(0, NUM_EXPERIMENTS):\n",
    "    result = dask.delayed(run_simulation)(\n",
    "        num_users=10000, \n",
    "        baseline_conversion=0.1, \n",
    "        relative_uplift=0.2, \n",
    "        is_aa_test=False,\n",
    "        randomization_unit='session',\n",
    "        persist_treatment_effect=False,\n",
    "        geometric_p=0.5,\n",
    "    )\n",
    "    delayed_results.append(result)\n",
    "\n",
    "results = dask.compute(*delayed_results)\n",
    "\n",
    "num_positive_effects_detected = np.sum([res['session_conversion']['test_gt_control'] > 0.95 for res in results])\n",
    "avg_effect_size = np.mean([res['session_conversion']['mean_diff'] for res in results if res['session_conversion']['test_gt_control'] > 0.95])\n",
    "\n",
    "print(\n",
    "    f\"Detected positive effect {num_positive_effects_detected/len(results):0.1%} of the time\\n\"\n",
    "    f\"Average effect size detected: {avg_effect_size*100*100:0.0f} bps\\n\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-24T21:35:54.095657Z",
     "iopub.status.busy": "2021-01-24T21:35:54.095197Z",
     "iopub.status.idle": "2021-01-24T21:38:36.389203Z",
     "shell.execute_reply": "2021-01-24T21:38:36.386747Z",
     "shell.execute_reply.started": "2021-01-24T21:35:54.095622Z"
    }
   },
   "source": [
    "## A/B Test - Session Level Randomization with Effect Persistence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T15:31:02.359860Z",
     "iopub.status.busy": "2021-01-31T15:31:02.359549Z",
     "iopub.status.idle": "2021-01-31T15:32:02.531898Z",
     "shell.execute_reply": "2021-01-31T15:32:02.527030Z",
     "shell.execute_reply.started": "2021-01-31T15:31:02.359812Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Detected positive effect 90.0% of the time\n",
      "Average effect size detected: 137 bps\n",
      "\n"
     ]
    }
   ],
   "source": [
    "delayed_results = []\n",
    "\n",
    "NUM_EXPERIMENTS = 250\n",
    "for iteration in range(0, NUM_EXPERIMENTS):\n",
    "    result = dask.delayed(run_simulation)(\n",
    "        num_users=10000, \n",
    "        baseline_conversion=0.1, \n",
    "        relative_uplift=0.2, \n",
    "        is_aa_test=False,\n",
    "        randomization_unit='session',\n",
    "        persist_treatment_effect=True,\n",
    "        geometric_p=0.5,\n",
    "    )\n",
    "    delayed_results.append(result)\n",
    "\n",
    "results = dask.compute(*delayed_results)\n",
    "\n",
    "num_positive_effects_detected = np.sum([res['session_conversion']['test_gt_control'] > 0.95 for res in results])\n",
    "avg_effect_size = np.mean([res['session_conversion']['mean_diff'] for res in results if res['session_conversion']['test_gt_control'] > 0.95])\n",
    "\n",
    "print(\n",
    "    f\"Detected positive effect {num_positive_effects_detected/len(results):0.1%} of the time\\n\"\n",
    "    f\"Average effect size detected: {avg_effect_size*100*100:0.0f} bps\\n\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-24T21:35:54.095657Z",
     "iopub.status.busy": "2021-01-24T21:35:54.095197Z",
     "iopub.status.idle": "2021-01-24T21:38:36.389203Z",
     "shell.execute_reply": "2021-01-24T21:38:36.386747Z",
     "shell.execute_reply.started": "2021-01-24T21:35:54.095622Z"
    }
   },
   "source": [
    "## A/B Test - User Level Randomization with Effect Persistence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-31T15:33:29.428042Z",
     "iopub.status.busy": "2021-01-31T15:33:29.427362Z",
     "iopub.status.idle": "2021-01-31T15:34:15.901476Z",
     "shell.execute_reply": "2021-01-31T15:34:15.880996Z",
     "shell.execute_reply.started": "2021-01-31T15:33:29.428008Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Detected positive effect 100.0% of the time\n",
      "Average effect size detected: 199 bps\n",
      "\n"
     ]
    }
   ],
   "source": [
    "delayed_results = []\n",
    "\n",
    "NUM_EXPERIMENTS = 250\n",
    "for iteration in range(0, NUM_EXPERIMENTS):\n",
    "    result = dask.delayed(run_simulation)(\n",
    "        num_users=10000, \n",
    "        baseline_conversion=0.1, \n",
    "        relative_uplift=0.2, \n",
    "        is_aa_test=False,\n",
    "        randomization_unit='user',\n",
    "        persist_treatment_effect=True,\n",
    "        geometric_p=0.5,\n",
    "    )\n",
    "    delayed_results.append(result)\n",
    "\n",
    "results = dask.compute(*delayed_results)\n",
    "\n",
    "num_positive_effects_detected = np.sum([res['session_conversion']['test_gt_control'] > 0.95 for res in results])\n",
    "avg_effect_size = np.mean([res['session_conversion']['mean_diff'] for res in results if res['session_conversion']['test_gt_control'] > 0.95])\n",
    "\n",
    "print(\n",
    "    f\"Detected positive effect {num_positive_effects_detected/len(results):0.1%} of the time\\n\"\n",
    "    f\"Average effect size detected: {avg_effect_size*100*100:0.0f} bps\\n\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}