733 lines
66 KiB
Plaintext
733 lines
66 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "ac27651a",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Exercise sheet\n",
|
|
"\n",
|
|
"Some general remarks about the exercises:\n",
|
|
"* For your convenience functions from the lecture are included below. Feel free to reuse them without copying to the exercise solution box.\n",
|
|
"* For each part of the exercise a solution box has been added, but you may insert additional boxes. Do not hesitate to add Markdown boxes for textual or LaTeX answers (via `Cell > Cell Type > Markdown`). But make sure to replace any part that says `YOUR CODE HERE` or `YOUR ANSWER HERE` and remove the `raise NotImplementedError()`.\n",
|
|
"* Please make your code readable by humans (and not just by the Python interpreter): choose informative function and variable names and use consistent formatting. Feel free to check the [PEP 8 Style Guide for Python](https://www.python.org/dev/peps/pep-0008/) for the widely adopted coding conventions or [this guide for explanation](https://realpython.com/python-pep8/).\n",
|
|
"* Make sure that the full notebook runs without errors before submitting your work. This you can do by selecting `Kernel > Restart & Run All` in the jupyter menu.\n",
|
|
"* For some exercises test cases have been provided in a separate cell in the form of `assert` statements. When run, a successful test will give no output, whereas a failed test will display an error message.\n",
|
|
"* Each sheet has 100 points worth of exercises. Note that only the grades of sheets number 2, 4, 6, 8 count towards the course examination. Submitting sheets 1, 3, 5, 7 & 9 is voluntary and their grades are just for feedback.\n",
|
|
"\n",
|
|
"Please fill in your name here:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"id": "d8b13e80",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"NAME = \"Kees van Kempen\"\n",
|
|
"NAMES_OF_COLLABORATORS = \"\""
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "e21d4bc9",
|
|
"metadata": {},
|
|
"source": [
|
|
"---"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "97ec4790",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "ba0b514985ed8d79af63814bda6ee328",
|
|
"grade": false,
|
|
"grade_id": "cell-70094cd1c9529a28",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"**Exercise sheet 3**\n",
|
|
"\n",
|
|
"Code from the lectures:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"id": "948639fc",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "1e25d190936687e34fe3ac01d4dfce44",
|
|
"grade": false,
|
|
"grade_id": "cell-9a6bf70e64f84715",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"import matplotlib.pylab as plt\n",
|
|
"\n",
|
|
"rng = np.random.default_rng()\n",
|
|
"%matplotlib inline\n",
|
|
"\n",
|
|
"def sample_acceptance_rejection(sample_z,accept_probability):\n",
|
|
" while True:\n",
|
|
" x = sample_z()\n",
|
|
" if rng.random() < accept_probability(x):\n",
|
|
" return x \n",
|
|
" \n",
|
|
"def estimate_expectation(sampler,n):\n",
|
|
" '''Compute beste estimate of mean and 1-sigma error with n samples.'''\n",
|
|
" samples = [sampler() for _ in range(n)]\n",
|
|
" return np.mean(samples), np.std(samples)/np.sqrt(n-1)\n",
|
|
"\n",
|
|
"def estimate_expectation_one_pass(sampler,n):\n",
|
|
" sample_mean = sample_square_dev = 0.0\n",
|
|
" for k in range(1,n+1):\n",
|
|
" delta = sampler() - sample_mean\n",
|
|
" sample_mean += delta / k\n",
|
|
" sample_square_dev += (k-1)*delta*delta/k \n",
|
|
" return sample_mean, np.sqrt(sample_square_dev / (n*(n-1)))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "bce418c0",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "90c9636c93064db3b008a2e9ec210117",
|
|
"grade": false,
|
|
"grade_id": "cell-8763c7e3277ef0ce",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"## Acceptance-rejection sampling\n",
|
|
"\n",
|
|
"**(35 points)**\n",
|
|
"\n",
|
|
"The goal of this exercise is to develop a fast sampling algorithm of the discrete random variable $X$ with probability mass function $$p_X(k) = \\frac{6}{\\pi^2} k^{-2}, \\qquad k=1,2,\\ldots$$\n",
|
|
"\n",
|
|
"__(a)__ Let $Z$ be the discrete random variable with $p_Z(k) = \\frac{1}{k} - \\frac{1}{k+1}$ for $k=1,2,\\ldots$. Write a function to compute the inverse CDF $F_Z^{-1}(u)$, such that you can use the inversion method to sample $Z$ efficiently. **(15 pts)**"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"id": "564c7f00",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "da96fc27e757b00d6e55befc89f5afda",
|
|
"grade": false,
|
|
"grade_id": "cell-b0b3fcf6d9bf35a0",
|
|
"locked": false,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def f_inverse_Z(u):\n",
|
|
" '''Compute the inverse CDF of Z, i.e. F_Z^{-1}(u) for 0 <= u <= 1.'''\n",
|
|
" assert 0 <= u and u <= 1\n",
|
|
" \n",
|
|
" # we need to round to return a natural number (excluding zero, yes).\n",
|
|
" return np.ceil(u/(1 - u)).astype(int)\n",
|
|
"\n",
|
|
"def random_Z():\n",
|
|
" return int(f_inverse_Z(rng.random())) # make sure to return an integer"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"id": "2ca9d9ca",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "3eeabc3f6572ed217a002adefda34876",
|
|
"grade": true,
|
|
"grade_id": "cell-a28cb07f90510b94",
|
|
"locked": true,
|
|
"points": 15,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"assert f_inverse_Z(0.2)==1\n",
|
|
"assert f_inverse_Z(0.51)==2\n",
|
|
"assert f_inverse_Z(0.76)==4\n",
|
|
"assert f_inverse_Z(0.991)==111"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "9a1b3dcd",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "78f978d46d5f742623cf174e88323ec9",
|
|
"grade": false,
|
|
"grade_id": "cell-576d3991eef6bf11",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"__(b)__ Implement a sampler for $X$ using acceptance-rejection based on the sampler of $Z$. For this you need to first determine a $c$ such that $p_X(k) \\leq c\\,p_Z(k)$ for all $k=1,2,\\ldots$, and then consider an acceptance probability $p_X(k) / (c p_Z(k))$. Verify the validity of your sampler numerically (e.g. for $k=1,\\ldots,10$). **(20 pts)**"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 5,
|
|
"id": "8e201507",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "862ecc31c9ee75e459dafe5262ddc776",
|
|
"grade": false,
|
|
"grade_id": "cell-20fd55e7a4c4c531",
|
|
"locked": false,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEICAYAAABI7RO5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAdWUlEQVR4nO3de3RU9dn28e9NCCcPoJxEUDmUBQICQgABJSAq1aJARNFlEaoQUUR43rdYxMqjVmsfYVkFBaRUBcRKAVEQWistKVpRQYwHVEAOlrwoIJUIIWAC9/tHknkIJCEbJ9mTyfVZKyuZPftw7z2Zuea3D79t7o6IiEhpVQm7ABERqVgUHCIiEoiCQ0REAlFwiIhIIAoOEREJpGrYBZSHevXqedOmTcMuQ0Skwvjggw++dff6RT0X18FhZtcC1/7kJz9h3bp1YZcjIlJhmNlXxT0X17uq3H2Zu6fWrl077FJEROJGXAeHiIhEn4JDREQCietjHFL+cnJyyMjI4NChQ2GXIiKlUKNGDZo0aUJiYmKpp1FwSFRlZGRwxhln0LRpU8ws7HJEpATuzt69e8nIyKBZs2alni6ud1WZ2bVmNiszMzPsUiqNQ4cOUbduXYWGSAVgZtStWzfwHoK4Dg6dVRUOhYZIxXEq79e4Dg4REYk+HeM4iVc2fh12CQCktGoUdgmnJNrbL+h2ePDBBzn99NP55S9/GbUarrnmGl566SUAXnrpJe66666ozHf8+PGsWLGCa665hsmTJ0dlnj/Gb3/7WyZOnBh2GSVat24dc+fOZerUqYGnffLJJ0lNTaVWrVrA/76uderUiXKVZSctLY0pU6bw+uuvl+ty1eIQCWjFihXUqVOHffv2MX369KjN99lnn2X9+vUxERqQFxyxIDc3t9jnkpKSTik0IC84Dh48GHlc8LrKyanFEcCkodeX6/Ienre4XJcXLx599FHmzp3LeeedR/369encuTMAW7ZsYfTo0ezZs4datWrxhz/8gdatWzN8+HDOPPNM1q1bxzfffMPjjz/O4MGD+frrrxkyZAjff/89ubm5zJgxg8suu4ymTZuybt06JkyYwJYtW+jYsSNXXnkl33zzDYMHD2bAgAEA3HLLLQwZMoTrrrsuUpu7c++99/KXv/wFM+PXv/51ZJysrCy6devGfffdx5AhQyLTvP/++4wbN47s7Gxq1qzJ888/T6tWrThy5Ai/+tWveOONNzAzRo4cyZgxY1i7di1jx44lKyuL6tWr8/e//51atWoxYcIE0tLSOHz4MKNHj+aOO+4gLS2NSZMmUbduXTZu3EivXr2YPn06EydOJDs7m44dO9K2bVvmz5/PwIED2bFjB4cOHWLs2LGkpqYCcPrppzN27Fhef/11atasyWuvvUbDhg3ZtWsXo0aNYuvWrQDMmDGDHj168OKLLzJ16lR++OEHunXrxvTp00lISCj0Gr7wwgssX76cQ4cOkZWVxbJlyxgzZgyffPIJubm5PPjggwwYMKDQN+6srKwixylqO7k7O3fupE+fPtSrV49Vq1ZFXtd69erxxBNP8NxzzwEwYsQIxo0bx/bt27n66qu59NJLeeedd2jcuDGvvfYaNWvWLFT7woULeeihh0hISKB27dqsXr2a7du3M3ToULKysgB4+umn6dGjB2lpafz3f/83DRs2JD09nZSUFC666CKeeuopsrOzefXVV2nRogXDhw+nRo0abNiwgV27dvHEE0/Qv3//Qsstbv03bNjAL37xC3744QeOHj3K4sWLadmy5Y97k7l73P907tzZT9XiL3ZGftp26V6uP8cuu6L47LPPCj0+dh2i8XMy69at83bt2nlWVpZnZmZ6ixYtfPLkye7ufvnll/umTZvc3f3dd9/1Pn36uLv7sGHDfPDgwX7kyBHfsGGDt2jRwt3dp0yZ4o888oi7u+fm5vr333/v7u4XXHCB79mzx7dt2+Zt27aNLDstLc0HDBjg7u779u3zpk2bek5OTqH6Fi1a5FdccYXn5ub6N9984+edd57v3Jm3XqeddlqR65SZmRmZz5tvvukpKSnu7j59+nRPSUmJPLd3714/fPiwN2vWzN9///1C0z777LP+m9/8xt3dDx065J07d/atW7f6qlWrvHr16r5lyxbPzc31K664whcuXFhkPXv37nV394MHD3rbtm3922+/dXd3wJcuXeru7uPHj48s58Ybb/Tf//73ke23b98+/+yzz7x///7+ww8/uLv7nXfe6XPmzDlhnZ9//nlv3LhxZJn33Xefz5s3z93dv/vuO2/ZsqUfOHDAV61a5T/72c9KHKeo7XTs61ig4HHB/9CBAwd8//793qZNG1+/fr1v27bNExIS/MMPP3R39xtuuCGyvGO1a9fOMzIyInW4u2dlZXl2dra7u2/atMkLPpNWrVrltWvX9p07d/qhQ4f83HPP9UmTJrm7+5NPPuljx45197z/0X79+vmRI0d806ZN3rhxY8/Ozi7V+t99993+4osvurv74cOH/eDBgyfUfPz71t0dWOfFfKbGdYvj2E4Oo0EtgNj31ltvMWjQoMh+64Jv+wcOHOCdd97hhhtuiIx7+PDhyN8DBw6kSpUqtGnThl27dgHQpUsXbrvtNnJychg4cCAdO3YscdnJycmMHj2a3bt388orr3D99ddTtWrht9jbb7/NzTffTEJCAg0bNiQ5OZm1a9cWapUcLzMzk2HDhrF582bMjJycHABWrlzJqFGjIss4++yz+eSTT2jUqBFdunQB4MwzzwTgb3/7Gx9//DGLFi2KzHPz5s1Uq1aNrl270rx5cwBuvvlm3n77bQYPHnxCHVOnTmXJkiUA7Nixg82bN1O3bl2qVasW+fbbuXNn3nzzTQD+8Y9/MHfuXIDIt+958+bxwQcfROrLzs6mQYMGRa73lVdeydlnnx2pf+nSpUyZMgXIO+373//+d6HxixunqO1UkrfffptBgwZx2mmnAZCSksJbb73FddddR7NmzSL/B507d2b79u0nTN+zZ0+GDx/OjTfeSEpKCpB3Yezdd99Neno6CQkJbNq0KTJ+ly5daNQo79hdixYtuOqqqwC46KKLWLVqVWS8G2+8kSpVqtCyZUuaN2/OF198Uar17969O48++igZGRmkpKT8+NYGcb6ryt2XAcuSkpJGhl2LlJ+iTi88evQoderUIT09vchpqlevHvk778sW9OrVi9WrV7N8+XKGDh3K+PHjufXWW0tc9tChQ5k/fz4vv/xyZFfHsQrmHcQDDzxAnz59WLJkCdu3b6d3796ReR2/rkUNKxg+bdo0+vXrV2h4WlraCeMXNX1aWhorV65kzZo11KpVi969e0fO/U9MTIxMk5CQUOIxCXdn2LBhPPbYY4WGL1myhIceegiA2bNnA0Q+uAumW7x4Ma1atSo0XUHIlzROcdukpBqLc+z/SUJCAtnZ2SeMM3PmTN577z2WL19Ox44dSU9PZ9q0aTRs2JCPPvqIo0ePUqNGjSLnWaVKlcjjKlWqFNqWJ3udilv/Cy+8kG7durF8+XL69evH7Nmzufzyy0vaBCelg+MSV3r16sWSJUvIzs5m//79LFu2DMj75t2sWTMWLlwI5L3JPvrooxLn9dVXX9GgQQNGjhzJ7bffzvr16ws9f8YZZ7B///5Cw4YPH86TTz4JQNu2bYusb8GCBRw5coQ9e/awevVqunbtWmIdmZmZNG7cGMjb91/gqquuYubMmZEPl//85z+0bt2anTt3snbtWgD2799Pbm4u/fr1Y8aMGZHWyqZNmyL7299//322bdvG0aNHWbBgAZdeeimQFwgF42dmZnLWWWdRq1YtvvjiC959990Sawbo27cvM2bMAODIkSN8//339O3bl0WLFrF79+5IzV999RWDBg0iPT2d9PR0kpKSTphXv379mDZtWuRD/cMPPyz1OEVtJyj69YO81+jVV1/l4MGDZGVlsWTJEi677LKTrm+BLVu20K1bNx5++GHq1avHjh07yMzMpFGjRlSpUoV58+Zx5MiRUs+vwMKFCzl69Chbtmxh69atJwREceu/detWmjdvzj333MN1113Hxx9/HHjZx4vrFoeEr7xPI+7UqRNDhgyhY8eOXHDBBYXe8PPnz+fOO+/kkUceIScnh5tuuokOHToUO6+0tDQmT55MYmIip59+emS3S4G6devSs2dP2rVrx9VXX83kyZNp2LAhF154IQMHDixynoMGDWLNmjV06NABM+Pxxx/nnHPOKXGd7r33XoYNG8YTTzxR6JviiBEj2LRpE+3btycxMZGRI0dy9913s2DBAsaMGRM5mL5y5UpGjBjB9u3b6dSpE+5O/fr1efXVVwHo3r07EyZM4JNPPqFXr14MGjQIgNTUVNq3b0+nTp147rnnmDlzJu3bt6dVq1ZccsklJdYM8NRTT5Gamsof//hHEhISmDFjBt27d+eRRx7hqquu4ujRoyQmJvLMM89wwQUXlDivBx54gHHjxtG+fXvcnaZNm0ZOQS345l3cOMVtp9TUVK6++moaNWpUaJdQp06dGD58eCTQR4wYwcUXX1zkbqmijB8/ns2bN+Pu9O3blw4dOnDXXXdx/fXXs3DhQvr06VOoNVVarVq1Ijk5mV27djFz5sxCrZaS1n/BggW8+OKLJCYmcs455zBp0qTAyz6enUrTuaJJSkryU72Rk67jCObzzz/nwgsvDLuM0Bw8eJCLLrqI9evXUxF6LAjrOoBoWbx4MUuXLmXOnDlhl1Kmhg8fTv/+/Ys89hQNRb1vzewDdz+x+Yd2VYlEzcqVK2ndujVjxoypEKFR0S1dupT777+fO+64I+xSKh21OE5CLY5gKnuLQ6QiUotDRETKVFwHh7pVFxGJvrgODle36iIiURfXwSEiItGn4JC4c2yvtWlpaSd0BlfWXnjhBXbu3Bl5PGLECD777LPA8wmjdpHSUHBI3Il2d+dFKalbjeODY/bs2bRp06ZM6xEpTwoOiTvHdnc+fvx4Dhw4wODBg2ndujW33HJLpEuGDz74gOTkZDp37ky/fv34+uu8U6/T09O55JJLaN++PYMGDeK7774DoHfv3kycOJHk5GSeeuqpIqdftGgR69at45ZbbqFjx45kZ2fTu3dvCk4H/+tf/0qnTp3o0KEDffv2BfK6/OjRowcXX3wxPXr0YOPGjSFsNZHSU5cjEnd+97vf8emnn5Kenk5aWlrkngTnnnsuPXv25F//+hfdunVjzJgxvPbaa9SvX58FCxZw//3389xzz3Hrrbcybdo0kpOTmTRpEg899FCk/6l9+/bxz3/+k5ycHJKTk4uc/umnn2bKlCkn9Lm0Z88eRo4cyerVq2nWrFmkz6TWrVuzevVqqlatysqVK5k4cSKLF6snZoldCg6JqmHDhkW6NC8PaWlpJx2na9euNGnSBICOHTuyfft26tSpw6effsqVV14J5HXC16hRIzIzM9m3bx/JyclA3voc2xV7wQ2WNm7cWOT0JXn33Xfp1asXzZo1A/63e+/iuk0XiVUKDol7x3eFnZubi7vTtm1b1qxZU2jck13zU9A5XXHTl6S47r2L6zZdJFYpOCSq5syZE3qXI8V1l32sVq1asWfPHtasWUP37t3Jyclh06ZNtG3blrPOOou33nqLyy67jHnz5kVaH6Wdvrjld+/endGjR7Nt27bIrqqzzz672G7TRWJVXAdHtO8AKBXDsd2d16xZk4YNG54wTrVq1Vi0aBH33HMPmZmZ5ObmMm7cONq2bcucOXMYNWoUBw8epHnz5jz//POBph8+fDijRo2iZs2ahVok9evXZ9asWaSkpHD06FEaNGjAm2++WWy36SKxSp0cnoQ6OQxGnRyKVDzq5FBERMqUgkNERAJRcEjUVYbdnyLx4lTerwoOiaoaNWqwd+9ehYdIBeDu7N2794T7l59MXJ9VJeWvSZMmZGRksGfPnrBLEZFSqFGjRuQC2dJScEhUJSYmRq6MFpH4pF1VIiISiIJDREQCUXCIiEggCg4REQkkroPDzK41s1kn6/FURERKL66Dw92XuXtq7dq1wy5FRCRuxHVwiIhI9Ck4REQkEAWHiIgEouAQEZFAFBwiIhKIgkNERAJRcIiISCAKDhERCUTBISIigSg4REQkEAWHiIgEouAQEZFAFBwiIhKIgkNERAJRcIiISCAKDhERCSSug0N3ABQRib64Dg7dAVBEJPriOjhERCT6FBwiIhKIgkNERAJRcIiISCAKDhERCUTBISIigSg4REQkEAWHiIgEUjXsAqR0Xtn4ddglAJDSqlHYJYhIyNTiEBGRQNTiqIAmDb2+XJf38LzF5bo8EYltanGIiEgganFUQGoBiEiY1OIQEZFAFBwiIhKIgkNERAJRcIiISCAKDhERCUTBISIigSg4REQkEAWHiIgEouAQEZFAFBwiIhKIgkNERAJRcIiISCAKDhERCaTCBYeZNTezP5rZorBrERGpjMo1OMzsOTPbbWafHjf8p2a20cy+NLMJJc3D3be6++1lW6mIiBSnvO/H8QLwNDC3YICZJQDPAFcCGcBaM1sKJACPHTf9be6+u3xKFRGRopRrcLj7ajNretzgrsCX7r4VwMxeBga4+2NA/1NdlpmlAqkA559//qnORkREjhMLxzgaAzuOeZyRP6xIZlbXzGYCF5vZfcWN5+6z3D3J3ZPq168fvWpFRCq5WLh1rBUxzIsb2d33AqPKrhwRESlJLLQ4MoDzjnncBNgZUi0iInISsRAca4GWZtbMzKoBNwFLQ65JRESKETg4zOy0/DOhAjOzPwFrgFZmlmFmt7t7LnA38AbwOfBnd99wKvMvYnnXmtmszMzMaMxOREQoxTEOM6tCXivgFqALcBiobmZ7gBXALHffXJqFufvNxQxfkT+vqHL3ZcCypKSkkdGet4hIZVWaFscqoAVwH3COu5/n7g2Ay4B3gd+Z2c/LsEYREYkhpTmr6gp3zzl+oLv/B1gMLDazxKhXJiIiMemkLY6C0DCzR45/ruBYR1HBIiIi8SnIwfHGZhY5RmFmDYCV0S8penRwXEQk+oIExx1Aqpl1NbMuwD+AKWVTVnS4+zJ3T61du3bYpYiIxI3SnFU1F1gPfAiMBl4CcoGB7v5l2ZYnIiKxpjQtjjn5491GXmg0Bb4Dfm5mg8uuNBERiUUnbXG4+9+Bvxc8NrOqQBugA3AJoBsqVSKvbPw67BJIadUo7BJEKrXS7Koyd490Oph/pffH+T/zihpHRETiV6kuADSzMWZW6KYWZlbNzC43sznAsLIp78fRWVUiItFXmgsAf0re8Y0/mVkzYB9Qg7w79P0N+L27p5dVgT+GuhwpO5OGXl+uy3t43uJyXZ6IFK80xzgOAdOB6flXiNcDst19XxnXJiIiMSjQjZzcPcfMbnX3/ymrgqRiUAtApPI6lftxfGlmT+Yf4/i5mekTRESkEjmV4FhCXktlZ/5vXcshIlKJnEpw/An4C9ATuAo4K6oViYhITAscHO4+xN2Xu/tG4L+AZ6JfVnTodFwRkegrdXCY2dVm9p6ZbTSzP5tZd3ffBYwow/p+FHVyKCISfUFaHNOB/0NeNyOzgMlmdrO7Z5VJZSIiEpOCnI67y93/lf/3SjNbA7xH3jEPERGpJIK0OLab2SNmVi3/cQ6wvwxqEhGRGBYkOBxIAXaY2dvAl0CambUsk8pERCQmlXpXlbvfDGBmNYB25HWr3gGYbWbN3f28silRRERiSaAuRyDSd9W6/B8REalkTuUCwApD13GIiERfXAeHruMQEYm+uA4OERGJPgWHiIgEouAQEZFAFBwiIhKIgkNERAJRcIiISCAKDhERCSTwleMiYXtl49dhlwBASqtGYZcgEoq4bnHoynERkeiL6xaHuy8DliUlJY0MuxYpG5OGXl+uy3t43uJyXZ5ILIrrFoeIiERfXLc4JP6pBSBS/tTiEBGRQBQcIiISiIJDREQCUXCIiEggCg4REQlEwSEiIoEoOEREJBAFh4iIBBLXwaG+qkREoi+ug8Pdl7l7au3atcMuRUQkbsR1cIiISPQpOEREJBAFh4iIBKLgEBGRQBQcIiISiO7HIXKKdO9zqazU4hARkUDU4hCJAt37XCoTtThERCQQtThEokAtAKlM1OIQEZFAFBwiIhKIgkNERAJRcIiISCAKDhERCUTBISIigcR1cOgOgCIi0RfXwaE7AIqIRF9cB4eIiESfgkNERAJRcIiISCAKDhERCUSdHIpUcLqhlJQ3BYdIHNF9QaQ8aFeViIgEohaHSBxRC0DKg1ocIiISiIJDREQCUXCIiEggCg4REQlEwSEiIoEoOEREJBAFh4iIBKLgEBGRQBQcIiISiK4cF5GoUGeLlYdaHCIiEohaHCISdeqlN76pxSEiIoGoxSEiUacWQHxTi0NERAJRcIiISCAKDhERCUTBISIigVS44DCzgWb2BzN7zcyuCrseEZHKplyDw8yeM7PdZvbpccN/amYbzexLM5tQ0jzc/VV3HwkMB4aUYbkiIlKE8j4d9wXgaWBuwQAzSwCeAa4EMoC1ZrYUSAAeO27629x9d/7fv86fTkQkIha6Pon3bk/KNTjcfbWZNT1ucFfgS3ffCmBmLwMD3P0xoP/x8zAzA34H/MXd1xe3LDNLBVIBzj///OisgIiIxMQFgI2BHcc8zgC6lTD+GOAKoLaZ/cTdZxY1krvPAmYBJCUleZRqFZEKRF2flI1YCA4rYlixH/TuPhWYWnbliIhISWIhODKA84553ATYGVItIhJHKksLoLzFwum4a4GWZtbMzKoBNwFLQ65JRESKUd6n4/4JWAO0MrMMM7vd3XOBu4E3gM+BP7v7higt71ozm5WZmRmN2YmICOV/VtXNxQxfAawog+UtA5YlJSWNjPa8RUQqq1jYVSUiIhWIgkNERAKJhbOqRETiSixcvQ5ldwV7XLc4dHBcRCT64jo43H2Zu6fWrl077FJEROKGdlWJiJSheOz2JK5bHCIiEn1qcYiIlKF47PZELQ4REQkkroNDZ1WJiERfXAeHzqoSEYm+uA4OERGJPgWHiIgEouAQEZFAFBwiIhKIgkNERAIxdw+7hjJnZnuAr8KuIwbUA74Nu4gYoW1RmLZHYdoecIG71y/qiUoRHJLHzNa5e1LYdcQCbYvCtD0K0/YomXZViYhIIAoOEREJRMFRucwKu4AYom1RmLZHYdoeJdAxDhERCUQtDhERCUTBISIigSg44pyZnWdmq8zsczPbYGZjw64pFphZgpl9aGavh11LmMysjpktMrMv8v9HuoddU5jM7L/y3yefmtmfzKxG2DXFIgVH/MsF/q+7XwhcAow2szYh1xQLxgKfh11EDHgK+Ku7twY6UIm3iZk1Bu4Bkty9HZAA3BRuVbFJwRHn3P1rd1+f//d+8j4YGodbVbjMrAnwM2B22LWEyczOBHoBfwRw9x/cfV+oRYWvKlDTzKoCtYCdIdcTkxQclYiZNQUuBt4LuZSwPQncCxwNuY6wNQf2AM/n77abbWanhV1UWNz9/wFTgH8DXwOZ7v63cKuKTQqOSsLMTgcWA+Pc/fuw6wmLmfUHdrv7B2HXEgOqAp2AGe5+MZAFTAi3pPCY2VnAAKAZcC5wmpn9PNyqYpOCoxIws0TyQmO+u78Sdj0h6wlcZ2bbgZeBy83sxXBLCk0GkOHuBS3QReQFSWV1BbDN3fe4ew7wCtAj5JpikoIjzpmZkbcP+3N3fyLsesLm7ve5exN3b0regc9/uHul/Fbp7t8AO8ysVf6gvsBnIZYUtn8Dl5hZrfz3TV8q8ckCJakadgFS5noCQ4FPzCw9f9hEd18RXkkSQ8YA882sGrAV+EXI9YTG3d8zs0XAevLORvwQdT1SJHU5IiIigWhXlYiIBKLgEBGRQBQcIiISiIJDREQCUXCIiEggCg4REQlEwSEiIoEoOERCYGZXmNm8sOsQORUKDpFwdCDvymSRCkfBIRKODsCHZlbdzF4ws9/m948kEvPUV5VIODoAu4E3gNnuXll76JUKSH1ViZSz/G7uvwW+Au5w9zUhlyQSiHZViZS/NsBa8npgPRJyLSKBKThEyl8H4B3y7gfyvJk1DLkekUAUHCLlrwPwqbtvAn4F/Dl/95VIhaBjHCIiEohaHCIiEoiCQ0REAlFwiIhIIAoOEREJRMEhIiKBKDhERCQQBYeIiATy/wEBQQ9uKI9JhwAAAABJRU5ErkJggg==\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"def accept_probability_X(k):\n",
|
|
" '''Return the appropriate acceptance probability on the event Z=k.'''\n",
|
|
" \n",
|
|
" # The resulting c is found by doing algebraic things, and seeing that\n",
|
|
" # 6/pi^2*(1 + 1/k) <= c should hold for all allowed k. For k = 1, c is\n",
|
|
" # the largest, so: tada.\n",
|
|
" c = 12*np.pi**(-2)\n",
|
|
" \n",
|
|
" # But we do not need this variable, we only need the acceptance probability.\n",
|
|
" return .5*(1 + 1/k)\n",
|
|
" \n",
|
|
"def random_X():\n",
|
|
" return sample_acceptance_rejection(random_Z,accept_probability_X)\n",
|
|
"\n",
|
|
"# Verify numerically\n",
|
|
"samples = [random_X() for _ in range(1000000)]\n",
|
|
"k = np.array(range(1,11))\n",
|
|
"\n",
|
|
"# The theoretical densities to compare the densities from the\n",
|
|
"# acceptance-rejection densities to:\n",
|
|
"p_X = lambda k: 6/(np.pi*k)**2\n",
|
|
"\n",
|
|
"plt.figure()\n",
|
|
"plt.hist(samples, k - .5, density=True, color=\"lightblue\", label=\"density of acceptance-rejection samples\")\n",
|
|
"plt.scatter(k[:-1], [p_X(i) for i in k[:-1]], s=800, marker=\"_\", color=\"black\", label=\"theoretical\")\n",
|
|
"plt.yscale(\"log\")\n",
|
|
"plt.ylabel(\"$p_X(k)$\")\n",
|
|
"plt.xlabel(\"$k$\")\n",
|
|
"plt.legend()\n",
|
|
"plt.show()\n",
|
|
"# TODO: The histogram overshoots the theoretical value every time."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 6,
|
|
"id": "1b0e1d47",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "641adac2c6462be641ae2d45a4d80e7d",
|
|
"grade": true,
|
|
"grade_id": "cell-04bad69f7cc2ad18",
|
|
"locked": true,
|
|
"points": 20,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"from nose.tools import assert_almost_equal\n",
|
|
"assert min([random_X() for _ in range(10000)]) >= 1\n",
|
|
"assert_almost_equal([random_X() for _ in range(10000)].count(1),6079,delta=400)\n",
|
|
"assert_almost_equal([random_X() for _ in range(10000)].count(3),675,delta=75)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "d010b3c8",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "ec5e4c478f45388c302b82613e45a0c0",
|
|
"grade": false,
|
|
"grade_id": "cell-f831d5b44932cf2f",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"## Monte Carlo integration & Importance sampling\n",
|
|
"\n",
|
|
"**(30 Points)**\n",
|
|
"\n",
|
|
"Consider the integral \n",
|
|
"\n",
|
|
"$$\n",
|
|
"I = \\int_0^1 \\sin(\\pi x(1-x))\\mathrm{d}x = \\mathbb{E}[X], \\quad X=g(U), \\quad g(U)=\\sin(\\pi U(1-U)),\n",
|
|
"$$ \n",
|
|
"\n",
|
|
"where $U$ is a uniform random variable in $(0,1)$.\n",
|
|
"\n",
|
|
"__(a)__ Use Monte Carlo integration based on sampling $U$ to estimate $I$ with $1\\sigma$ error at most $0.001$. How many samples do you need? (It is not necessary to automate this: trial and error is sufficient.) **(10 pts)**"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 7,
|
|
"id": "bbf5c843",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "febbb705202776d2af26bd95c4c08625",
|
|
"grade": true,
|
|
"grade_id": "cell-47eab0537d21690e",
|
|
"locked": false,
|
|
"points": 10,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"\n",
|
|
"After trial and error, I found that 43402 trials were sufficient to arrive at an 1-sigma error of around 0.001,\n",
|
|
"with mean 0.487692 and standard deviation sigma 0.000993. Wolfram Alpha tells me it should be approximately I =\n",
|
|
"0.487595, so I am happy.\n",
|
|
" \n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"def sample_X():\n",
|
|
" x = rng.random()\n",
|
|
" return np.sin(np.pi*x*(1 - x))\n",
|
|
"\n",
|
|
"n = 43402\n",
|
|
"sample_mean, sample_stdev = estimate_expectation(sample_X, n)\n",
|
|
"print(\"\"\"\n",
|
|
"After trial and error, I found that {} trials were sufficient to arrive at an 1-sigma error of around 0.001,\n",
|
|
"with mean {:.6f} and standard deviation sigma {:.6f}. Wolfram Alpha tells me it should be approximately I =\n",
|
|
"{}, so I am happy.\n",
|
|
" \"\"\".format(n, sample_mean, sample_stdev, 0.487595))\n",
|
|
"# TODO: Maybe I did not answer the question. Reading is difficult."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "26b8ff78",
|
|
"metadata": {},
|
|
"source": [
|
|
"__(b)__ Choose a random variable $Z$ on $(0,1)$ whose density resembles the integrand of $I$ and which you know how to sample efficiently (by inversion method, acceptance-rejection, or a built-in Python function). Estimate $I$ again using importance sampling, i.e. $I = \\mathbb{E}[X']$ where $X' = g(Z) f_U(Z)/f_Z(Z)$, with an error of at most 0.001. How many samples did you need this time? **(20 pts)**"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "e8ab8be3",
|
|
"metadata": {},
|
|
"source": [
|
|
"After hours of thinking, I tried $$\\sqrt{1.2x}$$ on $(0, 1)$. It works for the first half, so let's consider $$f_Z(x) = \\sqrt{1.2\\left( \\frac{1}{2} - |\\frac{1}{2} - z| \\right)}$$ on $(0, 1)$.\n",
|
|
"It has cumulative density $$F_Z(x) = \\int_0^x\\sec^2(\\pi{}t)dt = \\left[ \\frac{\\tan (\\pi{}t)}{\\pi} \\right]_{t=0}^x = \\frac{\\tan (\\pi{}x)}{\\pi}.$$ This density gives us $$X' := g(Z)f_U(Z)/f_Z(Z) = \\sin(\\pi{}Z(1-Z))\\sin^2(\\pi{}Z)$$ for $Z$ on $(0, 1)$. $Z$ can be sampled using the inversion method as $$F_Z^{-1}(p) = \\frac{\\arctan(\\pi{}p)}{\\pi}.$$\n",
|
|
"\n",
|
|
"Unluckily, I found that the function is way too non-constant."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 138,
|
|
"id": "8b01dfc9",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABioklEQVR4nO3dd1gUxxvA8e/QkSKi2AtiQcTeu2KLLbHEWKIxRpNoTDTWGGOLLVFTrTEaSxI1ll+MGms0omLHjqJYEBU7IL3D/P44JJQDTil3wHye5x693dmddwe493Z3dkZIKVEURVEUQ2Ok7wAURVEURRuVoBRFURSDpBKUoiiKYpBUglIURVEMkkpQiqIoikFSCUpRFEUxSCpBKYqiKAZJJShFURTFIKkEpRgEIYSfEKLjK257VQjRLmcjyrJOZyHEBSFEmBBijA7l1wkh5ib9P0fjTbm/7LRjVvvOSS/bfkrhpBKUkqGkD7tYIUSJNMsvCiGkEMJRT6GlIqV0lVIezuNqPwMOSyltpJSLX2ZDXePVNdnk1PFrqy8X2/aV208bIcQMIUR4mld00u9p/xyIV9EDlaCUrNwBBr54I4SoDVjqLxyDUQm4qs8AhBAm+qw/m165/bQdt5RytpTS+sULKA1cBPYA27ITqKI/KkEpWfkdGJLi/bvAby/eCCEmCSH+TLmBEGKJEOJHbTsTQlQQQmwTQjwTQgQKIZamWF1PCHFZCBEihNgshLBIsZ2LEOKwECI46bLTGynWpfrmn1EdQoiyQog/k5bfyezSUhb1HQLcgKVJ39Sra9m+vhDifNIlrM1AymNJG+9kIcSDpLI+QogOQojfgYrA30l1fJZi28lCiMtAhBDCRMuZT2MhhLcQ4rkQYm2adpRCiKop3q8TQszNor6OWbVJirITM/oZZtZ+Ou471XFn8rOzBHYBEcCbUsq4jMoqBk5KqV7qpfUF+AEdAR/ABTAG7qP59isBR6AMmg8Cu6RtTICnQEMt+zMGLgE/AFZoPrRbpajrDFAWsAeuASOT1pkCt4AvADOgPRAGOKeMM7M60HwZOwfMSNqHE+ALvKYlzkzrSypzGHg/g3YzA+4C45L21ReIA+Zqidc5qU3LJr13BKqkLZfmZ3IRqABYatmfH3Alab09cPxFvUnrJVA1xft12uLS8jugS5tk+DPU0kbJ7fcS+0513Jm0/V7gBGCt778h9creS51BKbp4cRbVCbgOPHixQkr5CDgKvJW0qAsQIKU8p2U/TdB8eE2SUkZIKaOllMdSrF8spXwopQwC/gbqJS1vBlgD86WUsVLKQ2i+IQ8kvYzqaAw4SM2loFgppS+wChigZR8vU582zdB86P4opYyTUv4P8MygbAJgDtQUQphKKf2klLez2P9iKeV9KWVUBuuXJq0PAua9RNyZ0bVNMvoZ5tS+MzzupLOqTUApoKuUMlzHY1MMlEpQii5+B94GhpLi8l4KvwKDk/4/OKm8NhWAu1LK+AzWP07x/0g0H1qgSTj3pZSJKdbfBcq9RB2VgLJJl5CChRDBaL6xl9Kyj5epT5uywAMpZcq5bO5qKyilvAWMBb4EngohNgkhymax//svsf5uUjzZpWubZPQzzIl9Z3jcQggjNGeDzkBnKWWIDvUqBk4lKCVLUsq7aDpLdEP7DeftQB0hRC2gB7Ahg13dByq+ws39h0CFpA+hFyqS4kxOhzruA3eklHYpXjZSym7ZrE+bR0A5IYRIs71WUsqNUspW/HfpdMGLVRltkkX9FdLU+zDF+0igSIr3pXXcb3bbJDO67juz+H5CcybWUUoZkAMxKQZAJShFV8OB9lLKiLQrpJTRwP+AjcAZKeW9DPZxBs2H93whhJUQwkII0VKHuk+juc/1mRDCVGiey3kdzeUcXes4A4Qm3Wi3FEIYCyFqCSEaZ7M+bU4C8cCYpE4MfdBcekxHaJ4Hai+EMAeigSg0l/0AnqC5V/ayPhZClBdC2KM5S9ycYt1F4O2k4+8CtE2xLrP6stsmmcnWvoUQ3wNdgQ5Jl5yVAkIlKEUnUsrbUsqzmRT5FahNxpf3kFImoPngqQrcA/yBLJ9RkVLGAm+g+RAKAJYDQ6SU13WtI8XyemjOBgOAX4Ci2akvk3j7oLkk+jzpGDPq6mwOzE+q5zFQEk1SAfgamJZ0SXKiLnUn2Qj8g6YTiC8wN8W6T9G0QzAwCM3Z7wsZ1pfdNslMdvaddNY+Ds2Z4FWR+jmoJ2nOypR8RqS+TK4or0YIURFNB4rSUspQfcejKEr+p75dKNmW9C11PLBJJSdFUXJKfn4SXTEAQggrNPcu7qLpYq4oipIj1CU+RVEUxSCpS3yKoiiKQdLbJb4SJUpIR0dHfVWvKIqiGIhz584FSCkd0i7XW4JydHTk7NnMei0riqIohYEQQutIK+oSn6IoimKQVIJSFEVRDJJKUIqiKIpBUglKURRFMUgqQSmKoigGSSUoRVEUxSCpBKUoiqIYJDUWn6LksUP3DnHU/2jye2NhTJ9qfXAt4arHqBTF8KgEpSh5JCYhhm88v2Gzz2ZszWyxMLYAIDwunG23tjGh4QQGuQwi9US8ilJ4qQSlKHngbuhdJh6ZyPWg6wx1HcqYBmMwNTIFICQmhGnHprHAcwGejz2Z3XI2Rc3TzaOoKIWOugelKLlsn98++u/qz6OIRyxtv5QJjSYkJyeAouZFWdx+MZMaTeKo/1H67+rPlYAreoxYUQyDSlCKkktiEmKYe2ouk45MoopdFbb22ErbCm21lhVCMMR1CL92/ZVEmcg7e99hw7UNqOlwlMJMJShFyQX3Q+/zzp532OyzmXdrvsu6LusoY10my+3qONRh6+tbaVW2FfPPzGf84fGExqpJipXCSSUoRclhB+8epN+ufviH+7PYbTETG09MdUkvKy8u+U1oOAH3++70/7s/3oHeuRixohgmlaAUJYfEJcSx4MwCxh0eh6OtI1tf34pbRbdX2pcQgqG1hrK2y1piE2N5Z887bPHZoi75KYWKSlCKkgMehT9i6P6hrL+2nkEug/it62+Usy6X7f3WL1mfra9vpXHpxsw5NYfPPT4nMi4yByJWFMOnEpSiZNOxB8d4a9db3A6+zXdtv+PzJp9jaqz7Jb2s2FvYs7zjckbXH80+v30M3D2Q28G3c2z/imKoVIJSlFeUkJjAkgtLGHVwFKWKlGJzj810duycK3UZCSM+rPMhqzqtIjgmmIG7B7Lbd3eu1KUohkKnBCWE6CKE8BFC3BJCfK5lfVEhxN9CiEtCiKtCiPdyPlRFMRyBUYGMODiClZdX0qtqLzZ020Al20q5Xm+TMk3Y+vpWXOxd+Nzjc+acnENsQmyu16so+pBlghJCGAPLgK5ATWCgEKJmmmIfA95SyrpAO+A7IYRZDseqKAbhwtML9Pu7HxefXmR2i9nMbjkbCxOLPKu/ZJGSrH5tNe/Veo8tN7YwZO8QHoQ/yLP6FSWv6HIG1QS4JaX0lVLGApuAnmnKSMBGaAYRswaCgPgcjVRR9ExKye/evzNs3zDMTcxZ3209vav11kssJkYmjG84nkVui7gXeo9+f/dLNQCtohQEuiSocsD9FO/9k5altBRwAR4CXsCnUsrEtDsSQnwohDgrhDj77NmzVwxZUfJeRFwEE49MZKHnQtqUb8OmHpuoYV9D32HRvmJ7NvfYTBmrMnz878csubCEhMQEfYelKDlClwSlbWjltA9jvAZcBMoC9YClQgjbdBtJuVJK2UhK2cjBweElQ1UU/bgdfJuBuwdy8N5BxjUcx49uP2Jrlu7XW28q2FZgfbf19Krai5WXVzLq31E8j36u77AUJdt0SVD+QIUU78ujOVNK6T1gm9S4BdwB9P/1UlGyad8dTbfukJgQfun8C8NqDTPI6TAsTCyY3WI2M5vPxPOxJ/129cPrmZe+w1KUbNElQXkC1YQQlZM6PgwAdqYpcw/oACCEKAU4A745Gaii5KW4RM2oEJOOTsK5mDNbemyhcenG+g4rU0II+lbvy+9df8cII97d9y5bb2xVo08o+VaWCUpKGQ98AuwHrgFbpJRXhRAjhRAjk4rNAVoIIbyAf4HJUsqA3ApaUXJTQFQA7+9/P3lUiDWvraGUVSl9h6Uz1xKubO6xmSalmzD75GxmnJhBdHy0vsNSlJcm9PXtqlGjRvLs2bN6qVtRMnLh6QUmHJ5AeFw4M5vPpLtT99ypKCFFJ1chwMg456tITOCnSz/x8+WfcbF34Qe3H3Jk+CVFyWlCiHNSykbplqsEpSiaLuQbr2/kW89vKWtdlh/cfqB6seo5XQncPACHv4aH51OsEODcFdp9DmXq5mydwJH7R5jiMQUjIyMWtl5Ii3ItcrwORckOlaAUJQNR8VHMPjmbXb67aFe+HfNaz8vZXnpSwu1/wf1reHAW7CpCnQFgnPQse3QwXPgdokOgRg9oNwVK18q5+oF7ofcYe3gst57fYnT90bxf+32D7OyhFE4qQSmKFvfD7jPOfRw3nt9gVL1RfFjnQ4xEDg5R6XsE3L+C+6egaAVoMxHqvg0maQZaiQqG0yvg5DKICYWavTSJqmTOdYaNjIvky5NfsvfOXtpXaM+8VvOwNrPOsf0ryqtSCUpR0jjx4ASTjk5CIpnfej5tyrfJuZ3fPQnu88DPA2zKQpsJUH9I+sSUVtRzTZI69RPERkDtvppEVbxKjoQlpWT9tfV8d/Y7KtpW5Ee3H3Eq6pQj+1aUV6USlKIkkVKy+spqFp9fTNViVVnUbhEVbCtkvaEuHpyHQ3M1l/SsSkLrCdBwKJi+5Fh9EYFwYhGcWQXxMVBvILSdrLk8mAM8H3sy8chEYhJimNdqHh0qdsiR/SrKq1AJSlHQXOaadnwaB+4eoItjF2a1mEUR0yLZ3/ETb80Z0/VdYGkPrcZB4/fBLJv7Dn8Kx34Az9UgEzXJrs1EsCmd7ZAfRzxmrPtYrgZe5cM6H/JxvY9z9vKmouhIJSil0Lsfep8x7mPwDfFlbIOxDHUdmv2OAoG3Nb3yvP4H5rbQYjQ0GwnmNjkT9AshD8DjWzj/GxiZQpMPNEmwiH22dhuTEMPcU3PZfms7bcu35evWX2NjlsOxK0oWVIJSCrVjD47x2dHPMBJGLGyzkBZls9nVOuQBHFkAF9aDiTk0HalJTtlMGFkKugOH58PlzWBmDS0+geYfZyshSinZ7LOZBWcWUN6mPIvcFuFkp+5LKXlHJSilUJJSsvbqWn489yPVilVjkdsiytuUf/UdRgaBx3eae0MyERq9B60ngk0ejzTx9JrmXtf1XVCkuCaGRsNe/l5XCueenGP84fHEJMTwVauvaF+xfQ4GrCgZUwlKKXSi4qOYeXwme/32Zv9+U0w4nFoOxxdDXITmOaZ2n0Ox3J9FN1MPzsG/s8H3MNiWB7cpUHfgK49MkfK+1Kh6oxhRZ4S6L6XkOpWglELlQfgDxrqPxSfIh08bfPrqo5DHx8L5X+HIQoh4qnmQtv30HH0+KUf4HoGDX2pGqHCooYmxRnfNMEovKTo+mtknZ/O37990qNiBea3mYWVqlfMxK0oSlaCUQsPzsScTDk8gPjGehW0X0qpcq5ffSWIiXN0Gh+bAcz+o1Ao6fgkVDHhEcynh2k74dw4E3oTyTaDTLKj08vfbUj4vVbloZRa3X0wFmxzqiq8oaWSUoNS5u1KgbL6+mQ//+RA7Czs2dt/4asnptjusagd/Dtd0RBj0Pxi6y7CTE2jOlmr2hFGn4PXFEHIf1naFjQM096xealeCd2q+w08df+Jp5FMG7h7IqUencilwRdFOnUEpBUJcQhxfn/marTe20qZ8G+a3nv/y3aUfXYaDM+H2IShaEdpPg9pvgVE+/R4XG6kZPunYjxAbBvXeBrepYFv2pXbzonv+nZA7TGo8ibdrvK3G8VNylLrEpxRYQdFBjHMfx/mn53m/9vt8Uu8TjF+mk0DwPTg0T9N129IO2kzSPGRrYv7SscTGJxIUEUtAeAzBkXGERccRFh1PWEw8MfEJxMYnEhufSEJi6r87E2OBmbEx5qZGWJgYYW1hio2FCTYWJthbmVHcyhx7KzOMjV4hMST3PFwJwhiaj4KWn4JFUZ13EREXwRSPKbjfd6dPtT5MbToVM+Mshm1SFB2pBKUUSD5BPow5NIbA6EBmt5hNN6duum8cFQzHvodTKzTvm32kefjV0i7DTaLjErgXFMmdgAjuBETg/zySh8HRPAyO4lFINCFRcVlWa2osUiUaKSE+UaZLWmkJAcWtzChrZ0mZohaUtbOkon0RHEtY4VTCinJ2lpgYZ3K299xP0zXda6uma3rbzzXd5I1Ns4wZIFEmsvTCUlZ5raJ+yfr80O4HilsW12lbRcmMSlBKgfPv3X+ZcmwKNmY2LHZbjGsJV902jI+Fs2s0D9pGBWm6jLefBnb/dQKQUvIwJBov/2C8H4Xh8zgUn8dh3A2KJOWfjF0RU8oWtUxOGg425hS3NqOEtTl2lqbYpDgTsjA1xszYCKMMzoLiExKJTUgkKjaB8Jh4wqLjCY2O43lEHIERMQSEx/IsLDo5IT4IjiIyNiF5e1NjQRUHa5xL2+Bc2oaaZWypU94Oe6s0ZzoPL8A/0zUD2dpX0XSkqNFD5x5/++7sY/rx6RSzKMbi9oupYW9gPRqVfEclKKXAkFKyymsVSy4soU6JOvzo9iMORRx02RCu74YD0yHIFyq3gU5zoGw9ouMSuHQ/GE+/IM7fC+ayfzAB4bEAGAlwLG6Fc2kbqpeywcnBisolrHAsYYWthW5nH7lBSklAeCx+gZqzudvPwrnxOAyfx2E8DPlvivcK9pbUKW9Hw4rFaFLZHpcythgL4MZ+ODADAnygYgt4bR6Ua6BT3VcDr/LpoU8JjQ3lq1Zf0bFSx1w6SqUwUAlKKRCi46OZcWIGe+/spbtTd2a1mIW5sQ73ih5egP1T4e5xKOFMfMdZXDBvgsetQE7cCuCyfwixCYkAVC1pTd3ydtStUJQ65e2oUdoGC9Ocn5I9N4VExXH1YQhe/iFc9g/h4v1gHgRHAWBtbkLDSsVoVbUErarY4fxwO0aHv4KIZ1CnP3SYAUWzHm3jWeQzxrqP5XLAZUbXH80HtT9QnSeUV6ISlJLvPYt8xphDY7gaeJUxDcYwvNbwrD8QQx9qRlq49AeJlsU55zSKXyJbc+z2cyJiEzASULu8Hc2c7GlcyZ5GjsWwK1Iwb/4/DI7C0y8IT78gTvkGcetpOAAlrM3oVKUIw9lOlVu/atq0xWhoORbMM5/QMCYhhpknZrLbdzddK3dldovZWJi8+nBLSuGkEpSSr3kHejP60GjCYsOY33p+1uPExUbCiSUkHvsBmRDPdouezHzehXCKUM7OErcaDrSq6kBzp+IULaK/y3T69CgkiuO3AvG4+YzDPs8IiYqjknEAXxf9ixaR7iRalcKo4wzNDMCZdLV/Mb/WovOLqFOiDovaL6KEZYk8PBIlv1MJSsm3Dt49yBfHvqCoeVGWtl+Ks71zxoWl5PmZPzA99CXWMU/YndCE+fEDsS9Xnc6upengUhLnUjbqUlQa8QmJnL8XzL/Xn7D/ymPsgy4y3XQ99Y1uEWTrgtnr32BdrXWm+3jRaaWoeVGWtF+iOk8oOlMJSsl3Un0zd6jDIreMv5lHxMRz0uMfKpyZg3OsN1cSHdlY7CMqN+xM19qlKV8sByYlLCSklFx7FMaeyw+IubCZ96J/pawI4qy1G9HtZtKsft0Mu7NfC7zG6EOjCY0N5evWX6uZehWdqASl5CuxCbHMOjmLnbd30q1yN2a3nK21M8SVByHsPH4Bl6s/0FscJoiinKs6muqvjaCSg60eIi9YpJRcufOIwH++odnj9SAl6417EdF4NH2aVqOCffrE/yzyGZ+6f8qVgCuMazguZyaGVAo0laCUfON59HPGuo/l/NPzfFzvY0bUGZHqAy4mPoHdlx+x/vhNGj3ezBiT7ViIOJ7VGkbpHtMQLzFCgqK7uKC7BPz1OWXu7+GBLM7X8YOIqNKDoa2caF21RKrnu6Ljo5l2fBr7/fbTp1ofpjWdhqmODwQrhY9KUEq+4Bviyyf/fsKTiCfMbTWXrpW7Jq97FhbD76fusvH0XVwjPZlrsZ4KiQ+Iq9IJ067zoURVPUZeiPgdJ3bXJMwCrnIOV6bGvENcCReGtnDkzYblKWJmAmhGnlh2cRkrL6+kSekmfN/ue4qaqy8PSnoqQSkG7/Sj04w7PA5TI1MWt19MXYe6ANwNjGDlUV+2nvOndMIjFtltoX7USaR9FUSX+VC9s54jL4QSE+D8r8h/Z0NUCLstuvNF8OsYFynGuy0cebe5I8WSRrD4+/bfzDwxk3LW5VjWYRkVbSvqOXjF0KgEpRi0v27+xeyTs3Es6sjSDkspZ12OG0/CWPzvTfZ4PcLaKI4fyh3CLXATRkYm0HYSNBv1SgO6KjkoMgjc58HZNcSZ2bHe5j1m+9fHwtSUgU0qMrKtEyVtLTj35Bxj3ccCsMhtEQ1K6TZihVI4qASlGKREmciSC0v4xesXWpRtwbdtv+XRc1iUlJiszIz50vkevR8vxjj0PtTqC53nvPSUEUoue3QZ9kyC+6eIKtWAn4p8xDIfa0yMBIOaVmJkWyeiecrH/37Mg/AHzG45mx5OPfQdtWIgVIJSDE5MQgxTj01lv99++lbvy+CqY1l00Je/Lz+kiKkxYxua8W7Icsx8D4CDC3T7Bipn/iyOokdSwqVNmrEOIwIIqz2Eb+L6seFyKCZGgqEtHBnUogQzT33G2SdntXaAUQonlaAUg/I8+jljDo3h4rOLjKj1KU/vN2XjmfuYGAs+aF6WUaa7sTz1o2YqiHZToOkInaeFUPQsKhgOf62Zf6pIcZ41n8bX/nX569JDbMxN+LBtJe4Z/cpev130rNKTmc1nqh5+hZxKUIrB8AvxY9S/o3ga+ZQ2xcaw/3RJouMTGdC4ApOqPMDu8BTNaOOufTQjbKvLefnTo8uwezz4e0KlltxuMouvPOHf608pZWtO0wbncH+ynqalm/K92/fYmqnn1gorlaAUg3Dh6QVG/zua+ETgyTAePy3Fa66l+KJ1MSp5zoWr26B4Vej2LVRx03e4SnYlJsLF9ZppPWLCoMVozlb6gFn77uD1IIQqla8TaLkeR9tKLO+4nLLW6stIYaQSlKJ3+/32M8XjC0S8HUG+7+Jc3JEZ3ZxpEbxTM+J4fAy0maiZjlz1zitYIgI1SerierCrSGKXb/gz3JWF+30ISvDGttIGbC0s+anjcmoWr6nvaJU8phKUojdSSlZdXsuSiz+QEOmIScAwJnVqwNsVQzDePRYenAMnN+j+HRSvou9wldzkd1xz2e/ZdajZk4j281jiGcHq0yewqLAOU9Movm/3Le0qttV3pEoeUglK0YuExATG/PMlR59sJy60Nl1LjWNqxyqUOPsDnFwGRezhta+hdl+dpxxX8rn4WDixGI4s1Jwpd5jBzQpv8fnfJ7mWuAhji4eMcJ3EJ43f0XekSh7JVoISQnQBFgHGwC9SyvlayrQDfgRMgQApZaZfgVSCKvgeh4YyaMenPE08S5Go9izqPJ1miZdg11gIvgcNhkDHWZokpRQ+gbdh1zi4cwTKN0G+vojN90356uwXSMvr1LXuy+o3pmGez2YzVl5eRgkq41nI/tvQGFgGdAVqAgOFEDXTlLEDlgNvSCldgbdyImgl/9p64TqdNw3iScI5mti+x9GBX9Ds4hewvg8Ym8PQPfDGEpWcCrPiVWDIDuj9MwTeQvzchgGhWzjY/2fKGrXlUvj/aL12BOfuBug7UkVPskxQQBPglpTSV0oZC2wCeqYp8zawTUp5D0BK+TRnw1Tyi5DIOD784x++9PwIaebPhLpzWF2lPOYrmsGVbdDmMxh5DBxb6jtUxRAIAXUHwCeeUKsPHF1IyfWd2df6bbqWe5co89O8s/tDFuy/RFxCor6jVfKYLgmqHHA/xXv/pGUpVQeKCSEOCyHOCSGG5FSASv7hcfMZHZdt5ETkl1haRLOm5Ve8570etn0A9k4w0gPaTwVTC32HqhgaqxLQZyUM/hPiYxDrurIw5jFT6k3GxOo2v96ZzBs/7eXW0zB9R6rkIRMdymi7c532xpUJ0BDoAFgCJ4UQp6SUN1LtSIgPgQ8BKlZUIxoXFDHxCSzc58O6C/9gVWE9JSzsWFO+M05b3weZCF3mQ5MPwUjdS1CyULUjjDoJh+bA6Z9526cc5VuNYKzPb/jHLqT7T4FMe60tg5tWVEMkFQK6nEH5AxVSvC8PPNRSZp+UMkJKGQAcBeqm3ZGUcqWUspGUspGDg8OrxqwYEN9n4fRZfoJfL/+FVcV1VLUtw+ZoE5wOzoPyjTQfNs0+UslJ0Z25NXRdAMP/ATMr2uyZzq+WLhQtkkCRSj8xc99eRq4/R3BkrL4jVXKZLgnKE6gmhKgshDADBgA705TZAbQWQpgIIYoATYFrORuqYmj+POdPjyXH8E/Yh2W5TTSyLM2v189R6tFVeH0xvLMdijnqO0wlv6rQBEYchVbjqX11NxsfPaO0hRlFnX7h8L2jdFvkgadfkL6jVHJRlglKShkPfALsR5N0tkgprwohRgohRiaVuQbsAy4DZ9B0Rb+Se2Er+hQdl8Dnf15mwtYLlKx0gET7nXSSlqzwPolNpZYw6hQ0fFc916Rkn6kFdJwJH/xLRYvi/H7zMlWFCZYVfkNan2XAylOsOuqLvp7nVHKXelBXeSn3AiP5aMM5rj58Tt36+/CN9mBAeBSfh8Zi3HU+1B2oEpOSO+JjweNbwj2+Y2yZ0pw2M8JR9MfLuz5dXEuz8K062FqoUdHzo1d+DkpRXnC//pQeSzy49zyY5o034xvtwcfPg/nCti7GH5+Cem+r5KTkHhMzcPsC6/cPsjzOhtfCI/CTm2nfxIMD1x7zxpJj3HiievkVJCpBKVmSUrLM/RbDfvWkrH0itZyX4B12iRnPIxjZZh5i8P/UlBhK3ilbH7MPj7Cg+hD6h4bjGbabbo23Eh4bS+9lx9l/9bG+I1RyiLrEp2QqMjaez/53mV2XH9GnZhS+cj73ZDQLcKBTz3VQrJK+Q1QKMXnvNCv2fshyi0RamTjw9PlnnPNPYFzH6oxuXxUjI3VGnx+oS3zKS3sUEsVbK06y2+sRcxp4czVuJo8So/ipQk86DflXJSdF70TFpnz03gmmWtfkeNxTbCw/Y4RLID8cvMGoDeeJik3Qd4hKNqgEpWh15UEIvZYdJyAwgN9r/M4v4WuJMjZmTetvaNphHhipXx3FQJhZMeDNzSyoMZRLxomcj5nH6pruHPB+SP+VJ3kaGq3vCJVXpD5llHT+ufqYt1acpB43WGw/g0mJXpiZ2bDu9a24Vu2m7/AURauuzSayuPVC/MzN+TFmF7vKf0PEE196LTvO9ceh+g5PeQUqQSmprD1+h1HrzzDDegeDTObwcVEoYV2W33ptx6m4s77DU5RMta7SjZVd1hFoYc1oy0B+tp5Gx/gj9P3pJEdvPNN3eMpLUglKASAxUfL13mus2eXOP7ZfUcxoO5+WcsCpWHV+fX0zZazL6DtERdFJ/VINWNt9A3GW9gwvVYyBxstZYraUMeuO8NcFf32Hp7wElaAU4hISmbj1Eo88fueA5VQuFHnG5JIO1ClZn9Vd12FvoeZsUvIXZ3tnfu32O2ZWJRleviJFjT35x/ILftvyP1Yeva1GnsgnVIIq5CJj4xm11oMWV6az2GwZW8o78aVdEZqXbcGKTiuwMbPRd4iK8kocizryW9ffKG5TjhHlynHTzpz/mc8ieP8C5u26QmKiSlKGTiWoQiwkKo4ZKzbyxb0P6W18jJ/rdecb41A6VerE4vaLsTSx1HeIipItZazLsK7LOioVrcxoO1MOO7flM9PNuJ0ZwdxNh0hQScqgqQRVSAWERbNh8TS+ChxH6SKJLGozjKUhXrxR5Q0WtlmImbGZvkNUlBxR3LI4q19bjYu9CxPi/NjV5iOamN5mlM9Qlv+ykth4NVOvoVIJqhB6/OQx3ot6MSpqBSHlWvJ9qwGsvf8P/Z37M6flHEyMdJnHUlHyj6LmRVnZeSUNSjXgi/t72NF9OkbWJRn98DP+WfwRUVHqWSlDpBJUIfPk+kkSV7ShedwZ7jT8nMU167Dp9naGug5latOpGAn1K6EUTFamVizvsJyW5Voy68rP7O06mtsV3qRH6Cbu/tCeqMB7+g5RSUN9GhUWUvL88HKKbeqBkPHcfn0zK4oGs/32DkbVHcX4huPVFNpKgWdhYsFit8V0qNiB+ee+40jTVpxtuJAKMbeIXdaKqOsH9B2ikoJKUIVBTDgRf7xHscNTOE1tnr29n5/D9rDnzh7GNhjLR/U+UslJKTRMjU35pu03dHXsyg/nfsCzYiKnOvyPJ/E2mG96i9iD8yBRjeFnCFSCKuie+RC3oh0WN3awWLyN9bDNrLy/iIP3DjK58WSG1x6u7wgVJc+ZGpnydeuv6VmlJ8svLueq9UVu99zB9oRWmB1bSPzvfSFSTSevbypBFWRX/yJxpRvhz58wUkyjzbDZ/HxjJkf8jzC92XQG1xys7wgVRW+MjYyZ3XI2b1V/i1Veq7gmd2LR92emxr+PvONB4s9t4MF5fYdZqKkEVRAlxMP+qbB1KNcSyvMWCxj53mCWXfuCEw9OMKvFLPo599N3lIqid0bCiOnNpjPAeQBrr67FK2Y9jd8cR9/YGQSFxyDXvAbnftV3mIWW6k9c0EQEwNah4OfBdtNufBkziJ/ea8Dya1/g+diTOS3n0LNqT31HqSgGQwjBF02/wMTIhPXX1jPAOZ5Bvd+h058OrLdbhevfY+DhBei6UDPtvJJnVIIqSB5ehM2DkeFP+a7IWH4Ja8bKIbVZ6TOV80/P81Xrr+jh1EPfUSqKwRFC8FnjzzA1MmXt1bX0d5aM6zmI13dYs6LMbjqfWwtPvaHfb2BTWt/hFhrqEl9BcXkrrHmNRJnI50UX8nNIU5a87cqa29M4//Q8X7f6WiUnRcmEEIJxDccxrNYwNvts5p7YwOddXfjw0etsrDgL+dgLfm4L/mf1HWqhoRJUfpeYAAdmwLb3SSzbgE9tf2DLIwcW9K3BhrszuPD0AvNbz6ebk5poUFGyIoRgbIOxyUnqidkffNS2Ml/cqMZal1VgagFru8HFP/QdaqGgLvHlZ9Eh8Of7cPMfZKPhjA8dyN+XnzKrVzV2PpmdnJy6Vu6q70gVJd94kaQEgtVXVvNWdcnAJr2YfeY+Jp3WMsR/JmwfCU+uQMdZYKw+RnOLatn8KsgXNvaHIF9k9++Z/bgZ2y/7Ma6zI0dC5qvkpCjZIITg0wafArD6ymr6VTeiW+2uzDjwCJu+i+nt8BOcXArPrkPfNWBRVM8RF0zqEl9+5HccVnWAiGcwZAero91Ye9yPd1qU4XLc95x7co6vWn2lkpOiZMOLJDXUdShbbmymXJX9NHOyZ9K2a3hUmwSvLwLfw7C6MwTd0Xe4BZJKUPnNhfXwW08oUhze/5e/Q5yYu/sar9Wy57HFCjwfezK35Vy6O3XXd6SKku8JIRjfcDyDXQazyecPatU+QpWSVny0/jxXy/SGd/6CsMfwSwe4e1Lf4RY4KkHlF4mJcPBL2PExOLaE9w9wMtiOCVsu0dDRGkr+yulHp5jdcjavV3ld39EqSoHxogv62zXeZvONDbRqchprC2PeW+uJv10jeP9fsLCD396Ay1v0HW6BohJUfhAXDX8Oh2M/QMOhMOh/3AozZcTvZ6lQ3JySVbZy4tFxpjefTq+qvfQdraIUOEIIPm/yOf2q92PLzd/p1tqLqLgE3lvrSah1JXj/IFRoCts+gCPfgFQz9eYElaAMXUSg5pLe1W3QaTb0+JGgaMmwdZ6Ymkiq19rJsYdHmNJkCm9Vf0vf0SpKgSWEYGqzqfSq2outt9fQq9117gRE8PGG88Sb28HgbVBnALjPhZ2fQEKcvkPO91SCMmRBd2B1J80wK2+tg5afEpOQyIjfz/I4NJJGjQ5y7NG/TGw0kbdd3tZ3tIpS4BkJI75s/iXdKndj+91feL31TTxuBvDl31eRxqbQewW0/Vxzr3hDX4gO1XfI+ZpKUIbq4QVNcooKgnd3gmtvpJR8/qcXnn5BtG52lBNP9vNJvU941/VdfUerKIWGsZEx81rNo1OlThx8+gsdGvuy/tQ91h73AyHAbQr0XA53PGBddwh7ou+Q8y2VoAzRrX9hXQ8wsYBh/0DFZgAsc7/FXxf8adH4BGcC9zC81nBG1B2h52AVpfAxMTJhQesFtC7XGs/wVdR38WXubm/cfZ5qCtQfBG9vgcDbmi+aAbf0G3A+pRKUobm8BTb2g2KOMPwAOFQH4ID3E7795wa1a53GK/xvBtYYmPwgoaIoec/U2JTv231Po9KNuCNWU7GCL2P+uMDtZ+GaAtU6wtC/ITZCk6T8z+k34HxIJShDcmqFphdQxebw3h6wLQPAjSdhjN10gUqVPfFL2E6vqr34vMnnapp2RdEzCxMLlrRfgmsJV0Js1mJidZMPfj1LSFRSB4lyDWH4P2BhC7++Drfd9RtwPqMSlCGQEty/gn2ToUYPGPS/5KFTgiNj+eC3s5gVO0uQxZ90rtSZL5t/iZFQPzpFMQRWplYs77CcKkWdMCr9K/6R3oz54wIJiUldzYtXgWH7NVdFNvaDq9v1GW6+oj7l9C0xEfZMgiMLoP5geOtXzYjJQHxCIqP/uMCT+DMk2G+lZdmWzG89H2MjYz0HrShKSkXNi7Ki0wpKW5XEtvJveNy9zMJ91/8rYFMa3tsNZetrJhQ9t05foeYrOiUoIUQXIYSPEOKWEOLzTMo1FkIkCCH65lyIBVhCPOwYBZ6roMVoeGNpqpGRv/3nBicensC87CbqlazL9+2+x9TYVI8BK4qSkRKWJVjVeRVFLaywr7KOlSfPsMfr0X8FLItphkaq2gH+/hROLtNfsPlElglKCGEMLAO6AjWBgUKImhmUWwDsz+kgC6T4WM3oEJf+ALdp0GmOpotqkn1XHrHytDs2FTdQrVgVlnZYShHTInoMWFGUrJS1LsuqTquwNDPCzmktk7Yd5dbTsP8KmFnBgD/A5Q3Y/wUc/UZ/weYDupxBNQFuSSl9pZSxwCagp5Zyo4E/gac5GF/BFBcNW94B7+3QeR60nZQqOd16Gs7E7f9g4/grZaxLsqLTCmzNbPUXr6IoOnOyc2JFx58wNYvCuOwvvL/eg7DoFKNKmJhB37VQpz8cmgsHZ6mhkTKgS4IqB9xP8d4/aVkyIUQ5oDewIrMdCSE+FEKcFUKcffbs2cvGWjDERcGmgXBjH3T/Dlp8kmp1eEw8H2w8gCjzC3YWFqzs/DMlLEvoKVhFUV6FawlXFrn9iJF5AE8tf2LCVk9kyiRkbAK9VkCDd+HY97B/qkpSWuiSoLT1ZU7bkj8Ck6WUCZntSEq5UkrZSErZyMHBQccQC5C4KPhjoKarac9l0Pj9VKullEz433GeFlmChXkcKzv/TAWbCnoKVlGU7Ghetjnz23yNcZG7eIT8wM9Hb6YuYGSkmVOqyQg4tUxzyU8lqVR0mVHXH0j5KVkeeJimTCNgU9JzOSWAbkKIeCnl9pwIskCIjdScOfkegV7LoV76sfPWnLyBR9gCzCyD+anjSpztnfUQqKIoOaWLYxeCo4OZd3oeiy59TaNK39HI0f6/AkJA1wWaf08t1ySoLl+nuuRfmOlyBuUJVBNCVBZCmAEDgJ0pC0gpK0spHaWUjsD/gFEqOaXw4rKe7xHo9ZPW5HTpfiDfXZiOsYU/37RdSKPSjfQQqKIoOW1AjQEMrfkBJkXP8uGueQRFxKYuIAR0mQ9NP4LTP8G+KepMKkmWCUpKGQ98gqZ33jVgi5TyqhBipBBiZG4HmO/Fx8CmQZrk1HsF1BuYrkhIVCzDdk3B2Poa4xpMppNjRz0EqihKbhnfaDRuZXsQZ/MP72z5kcTENAlICM2ZU7OPNUnqwAyVpNDtEh9Syj3AnjTLtHaIkFIOzX5YBURCHGx9D27/q3nGqe6AdEWklLy9dTaxRU7yesV3GFZnkB4CVRQlNwkh+L7DHN7c9oTb4euZvK8c33QblLYQvDYPEmLgxGIwtQS3L/QTsIFQI0nkloR4zbh6Pruh27fQ4B2txT775xfuyR24WLdnXrtJeRykoih5xcTIhI1vLMHWqDJ7n3zL5sse6QsJAV2/0Ywqc2QBeHyf94EaEJWgckNiIuwcDVf/0jyA2+QDrcU2XTnA3kdLsZWu/N7zGzX4q6IUcFZmVmx6YxXGiXbMOzsJ72e+6QsZGcHri6FWX/h3lmYQ6UJKJaicJiX8Mw0ubdTMrNlyjNZil55489XZKRjFlWZjz2WYm5jlcaCKouhDRbuSfN1yEYlIhu75kKCooPSFjIw196xr9NAMIn1pc94HagBUgsppx77XPNPQZAS00z5s4eOIx7y//yMS4s2Z2/wHKhUrnsdBKoqiT91q1KF32WlEJgYy6O+PiEmISV/I2BTeXA2OrTVjdt4ofKPIqQSVk86uhX9nQ+23NN1GtVyyi4iLYMjuEUTFR9C91HTeqJVuWENFUQqBWZ17UClxOP5R3nx6cDKJMjF9IVMLGLARSrnCliFw92TeB6pHKkHllGt/w65xULWT5lkno/RNG58Yz5h/J/Ao0o/SMR8yt2tnPQSqKIohMDISrHnrfYye9+D4439ZfH6p9oIWtjDoTyhaHjb2hyfeeRuoHqkElRPun4E/39fMntnvN82puRbfeH7DmSfHSXjWk5V938bMRDW/ohRmpWwtWNBxDLHBjVh9ZRU7b+/UXtDaQTNVh6klbHgLQh9pL1fAqE/I7Aq8rflWY1MG3t4MZtqnxNh4bSMbr28kNrAVX7QejpODdR4HqiiKIepSqwxvlB1NfEQVZhyfydnHZ7UXtKsIg7ZA1HPY+BbEhGkvV4CoBJUdEYGwIWluxsF/gpX2UcePPzjOgjMLSAh3oYX9UAY1rZiHQSqKYuhmvlGH4hEfIOPs+dR9LPdD72svWKYu9PtVc5lv61DNYAAFmEpQryouGv4YAKEPNWdOxatoLeYb4svEIxMxTiiDxfMhLHyznnreSVGUVKzNTVjUrwXhd98lKjaBTw59QnhsuPbC1TpBj+/h1kHYM7FAD4mkEtSrkFIzZbP/Gej9M1RoorVYSEwIo/8dTXyCEc99B7OgT2McbMzzOFhFUfKDhpWK8XGrpoTcHYhfyF0+O/oZCYkZzGDUcCi0Gg/n1sHpn/MyzDylEtSrOLEYLm8Ct6ng2ktrkfjEeCYemciD8Ic8v/M2b9atTaeapfI2TkVR8pXRHapRw64+IqgXHg88WHR+UcaF208H5+6wfwrc+jfvgsxDKkG9LJ99cGAmuPaBNhmPnfeN5zecenQKy9B+lDStwYzX1fNOiqJkztTYiO/61SUyoCllRHvWXl3Ljls7tBc2MoI+K6FkTc2g1AE3tZfLx1SCehlPr8Gfw6FsPc2MuBncS/rr5l9svL6R6hbdeeRfhwV962Brob3ruaIoSko1StsytlM1bnh3wMm6HrNPzuZKwBXthc2tYeAfmkdbNvbX9PArQFSC0lV0iGZeJzMrzZPdGXQnv/zsMnNOzcHFriHnL7ZgUNOKtK5WCKe3VxTllX3Y2ol6FYpzx7s3xSyK86n7pwREBWgvbFcR+q+H4Lvw10jNYNUFhEpQupASto+C537w1jqwLau1WEBUAOPcx1HC0oFHN/tSvpg1X3RzydNQFUXJ/0ySLvXFxFjiEDmC0JhQJhyeQFxG3corNYfXvoIb++DYd3kbbC5SCUoXxxfB9V3QeQ5UaqG1SFxCHOMPjyc0NpTaJp9yP0Cw4M06WJnrNCekoihKKlUcrJnY2ZmT1yx4vdynnH96ngWeCzLeoMmHmik6Ds2D24fyLtBcpBJUVnyPaOZkce0NzUZlWGyh50IuPL3Ae86T+eu05O2mFWlRRfuDu4qiKLoY1qoydSvYsf1YKfpXe4fNPpsz7jQhBLyxGBxqwP+GQ3AGD/vmIypBZSbsMfxvGBSvCm8sybBTxC7fXWzy2cSgGu+wzcOBUrYWTOlaI4+DVRSloDE2EnzTtw5h0XE8uuNGk9JNmHNqDteDrmvfwMxKcz8qIQ62vpvvR5pQCSojiYnw1wiIjdAMAGtuo7XYzec3mX1yNg1KNsA4uDu3nobzVZ/a2Khee4qi5IDqpWwY3b4au72e0q3URIqaFWWc+zhCY0O1b1CiKvRcCg/Ogfu8vA02h6kbJBk5uQR8D0OPH6Ck9o4O4bHhjDs8DitTKz6sMZN3V12nT4NyuDmXzNtYk8TFxeHv7090dLRe6leUgs7CwoLy5ctjapq3X0A/aleFvVces2D3A75/ZwGjD3/AVI+pLGq/CCOh5TzDtRfcHgLHfgQnN3Bqm6fx5hQh9TSOU6NGjeTZsxmM2qtvD87D6k7g3BX6/a710p6UkvGHx+N+351VnX5hzp/RPAyO4uD4ttgV0c/07Xfu3MHGxobixYur8f4UJYdJKQkMDCQsLIzKlSvnef1XHoTwxtJj9G9ckVoul5l/Zj5j6o/hgzofaN8gNgJWtoPoUPjoBFgZ7szdQohzUspGaZerS3xpxYRpHsa1LgWvL87wvtPG6xs5eO8g4xqOw+t2cS77hzDjdVe9JSeA6OholZwUJZcIIShevLjerlDUKleU4a0q88eZe1Q170LXyl1ZenFpxtNzmFlB3zUQFQQ7Ps6Xg8qqBJXWvima5536rIIi9lqLXAm4wrdnv6Vd+Xa0L9OXb//xwc3ZgdfrlMnbWLVQyUlRco++/77GdapOOTtLvvjLi88bT6OCTQUmH51MUHSQ9g1K14ZOc+DGXji7Om+DzQEqQaV08yBc+B1afgqOLbUWCY0NZeKRiThYOjCn5Rxm7vRGSpjTq5bef3kVRSnYipiZMK93LW4/i+C344/5ru13BMcE84XHFyTKDEaQaDpCcx/qnxnw/G7eBpxNKkG9EB0Cf4/RPEPQborWIlJKZh6fyZOIJyxss5BjNyI5dP0pEzpXp3wx7UMfKYqi5KR2ziXpWa8sy91vYxxflslNJnP84XHWXFmjfQMhkh6TMYKdn+SrS30qQb2wfyqEPYKey8FE+5xNm3w2cfDeQT5t8CmVbVz5cqc3dcoX5b2WeX/DND+Kioqibdu2JCRkMMcNEBsbS5s2bYiPj8/Rulu00D4CSEq6xPcy2w8bNoySJUtSq1atV9pfZtK2U261m2KYpveoiaWZMV9su0Lfan3p4tiFpReWcv7Jee0b2FXQjIRz5yiczSCRGSCVoCD1pb3yDbUXeX6Tbz2/pVW5VgxxHcK3+30Iiojhq961MTZSl/Z0sWbNGvr06YOxsXGGZczMzOjQoQObN2/O0bpPnDiRZRld4nuZ7YcOHcq+ffteaV9ZSdtOudVuimEqYW3OlK41OOMXxJ/nHzCz+UzKWJVhiseUjJ+PajgUnNrBgfxzqU8lqOhQzaW9Es7Q9nPtReKj+ezoZ9iY2TC35Vy8/ENZf/ouQ5o7Uqtc0TwO2PBdu3aNNm3aUKdOHb755huqVq0KwIYNG+jZs2dyOTc3Nw4cOADAtGnTGDNmDAC9evViw4YNr1R3REQE3bt3p27dutSqVSv5A9va2hoAPz8/XFxc+OCDD3B1daVz585ERUWli+/WrVs4ODjg6OhIvXr1sLe3p0qVKoSGhmYYd9rja9OmDfb22jvavCxtdaZtp+y0m5L/9GtUgQYV7fh673Xi4sxY0GYBTyKfMPfkXLQ+PvTiUh/AztH54lKfelD38HwIfQjDD4CphdYiP5z7gVvBt/ip40/YmdszdPtxSlibM75z9TwOVnez/r6K98MMvkm9opplbZn5umumZeLj4xk0aBCrV6+mfv36fPTRR9SqVYvY2Fh8fX1xdHT8L8ZZs5gxYwZPnz7lwoUL7Ny5E4BatWrh6emZbt+tW7cmLCws3fJvv/2Wjh07ArBv3z7Kli3L7t27AQgJCUlX/ubNm/zxxx+sWrWKfv368eeff9KvX79U8VWtWpVWrVoxfvx4WrduTbt27ViyZAm2trZa49Z2fDkpo7ZK2U4ZtZtSMBkZCeb1rk2PJcdYuP86X/epw6h6o1hyYQmtyrfijSpvpN/IriJ0mgW7J8CVP6F237wP/CUU7gT15CqcXgEN34UKjbUWOep/lI3XNzLYZTCtyrXi1xN+eD0IYfHA+moSQi22bdtG3bp1qV+/PgA1a9akZMmSBAQEYGdnl6psmzZtkFLy/fffc/jw4eRLY8bGxpiZmREWFoaNzX9DTHl4eGRZf+3atZk4cSKTJ0+mR48etG7dOl2ZypUrU69ePQAaNmyIn5+f1viuXr2afP/o+vXrODs7Zxj3kydP0m2flY4dO/L48eN0y+fNm5fqTCyjOoFU7ZRRuykFl0sZW95r4cgvx+7Qt2EFhtcazomHJ5h3ah71HOpR0bZi+o0avgcX1mvuu1frDBa2eR+4jgpvgpIS9kzS/HA6zNRaJDAqkOnHp+NczJlxDcfxNDSab/f70KpqCYN45ikzWZ3p5JbLly8nf/gDXLlyhS5dumBpaZnuAUcvLy8ePXpEiRIl0n2gxsTEYGGR+oxWlzOo6tWrc+7cOfbs2cOUKVPo3LkzM2bMSFXe3Py/TjDGxsZERUWliy8qKoro6GiKFSvG/fv3KV68OGZmZhnGre34snLw4EGdy2bUVmnbSVu7KQXb2E7V2e31iGnbr/D3Jy2Z33o+fXb24XOPz/mt62+YGKX5mDcyhu7fwaoOcGQBvGa44/UV3ntQXlvh7nHo+KXWB3KllMw6OYvw2HDmt56PmbEZX++9Tkx8IrN7uqpnnjJQvHhxbty4AcDFixdZv349devWpVixYiQkJCR/iD969IhBgwaxY8cOrKys2L9/f/I+AgMDcXBwSDfemYeHBxcvXkz3epGcAB4+fEiRIkUYPHgwEydO5Pz5DHo1pZE2Pm9vb1xcNGMwXrt2Lfn/GcWddvuclFGdadspo3ZTCjZrcxNmvl6Ta49C+e3kXUpblWZG8xl4BXix2iuDh3PLNdRcOTr1EzzxztuAX0LhTFDRIfDPNM0Pqf4QrUV23t6J+313xjQYQ9ViVTlzJ4i/LjzgwzZOODlY53HA+cc777zD2bNnady4MWvWrMHR0REnJycAOnfuzLFjx4iMjKRPnz589913uLi4MH36dL788svkfbi7u9OtW7dXqt/Ly4smTZpQr1495s2bx7Rp03Te9kV8kPrynqWlJefPn8fb2zvTuFNuDzBw4ECaN2+Oj48P5cuXZ/Xql3+SP7O2SttO2Wk3JX97zbU0bao78MOBGzwLi6GLYxe6OnZlxaUVGU/N0WGm5grSnomG22FCSqmXV8OGDaXe7J0i5cyiUvqf07r6YdhD2WxDMzlkzxAZnxAv4+IT5Gs/HJHNvzooI2Li8jbWl+Dt7a3vEGRYWFjy/xcuXCinTp2a/P78+fNy8ODBWe6jd+/e8vr167kSX2Z0jS+3tn9ZadtJX+1W2BjC35k2t56Gyapf7JYTtlyUUkoZHB0s3Ta7yV7be8mY+BjtG51dK+VMWykvbcm7QLUAzkoteaLwnUE9vwtnVkL9wVCuQbrViTKR6SemkyATmNtqLsZGxqw/dZfrj8OY1qMmRcwK7207Xfzwww+4urpSr149/Pz8mD59evK6+vXr4+bmluWDur169UrukJCXdIkvN7d/GWnbSZ/tphiGKg7WDG/lxP/O+XPu7nOKmhflyxZfciv4FssuLtO+Uf0hUKYuHJoN8bF5G7AutGWttC+gC+AD3AI+17J+EHA56XUCqJvVPvV2BrVtpJRzSkoZ7K919cZrG2WtdbXkFh/NN4pnYdGy1sx9ctCqUzIxMTEvI31phvrNTlEKEkP+OwuPjpNN5x2U3RYdlfEJms+rmcdnyjq/1pEXnlzQvtHNg5qzqFM/512gafCqZ1BCCGNgGdAVqAkMFELUTFPsDtBWSlkHmAOszInkmeOeeMOlP6DJB1C0XLrVD8Mf8sO5H2hZtiV9q2meD1i47zpRsQl8+UZN1TFCURSDZmVuwtTuLlx9GMofZ+4BMKnxJEoXKc2MEzOISYhJv1GV9uDYGo4uhJjwPI44c7pc4msC3JJS+kopY4FNQKqHNKSUJ6SUz5PengLK52yYOeTQXM3U7a3Gp1slpWT2ydkAzGg+AyEEF+8Hs+WsP8NaVaZqSfVciaIohq9HnTI0dyrON/t9eB4Ri5WpFTObz+ROyB1WXtZy7iCEpsNExDNNrz4DokuCKgfcT/HeP2lZRoYDe7MTVK647wk+u6HlGK3dynf57uL4w+OMbTCWstZlNd3M/75KCWtzRrevqoeAFUVRXp4Qgplv1CQsOo4fD2oe+WhRrgVvVHmDNV5r8AnySb9RhcZQowecWAyRGcwtpQe6JCht17W09kkUQrihSVCTM1j/oRDirBDi7LNnz3SPMrukhINfgpUDNP0o3erAqEAWeC6gnkM9BtQYAMCOiw+5cC+Yz7o4Y6NGjFAUJR+pUdqWQU0rsf70PW480TzcPqnRJGzNbZlxYgbxiVpGvW8/TTOj+LHv8zjajOmSoPyBCinelwcepi0khKgD/AL0lFIGatuRlHKllLKRlLKRg4PDq8T7au4cgbvHoM1nYJ7+Gab5Z+YTGRfJrBazMBJGRMTE8/Xea9QuV5S+DQzzaqWiKEpmxneqjrW5CXN2eSOlxM7CjilNp+Ad6M167/XpNyjpAnUHwumVEPYk7wPWQpcE5QlUE0JUFkKYAQOAnSkLCCEqAtuAd6SUN3I+zGw6sRSsSkKD9A/lHvU/yj6/fYyoMwInO80DpSuO3OZJaAxfvlETIzWVhqIo+VAxKzPGdayGx80ADl57CsBrlV7DrYIbSy8u5UH4g/QbtZkICbHg+UseR6tdlglKShkPfALsB64BW6SUV4UQI4UQI5OKzQCKA8uFEBeFEGdzLeKX9fQ63DoATT5MN1p5TEIMX5/+mspFKzOs1jAA7gdF8vNRX3rWK0vDSjkzVYKiKIo+DGpWiWolrZm725uY+ASEEHzR9AuMhBELzixIv0HxKuDcVZOgYiPzPuA0dHpQV0q5R0pZXUpZRUo5L2nZCinliqT/vy+lLCalrJf0apSbQb+UU8vAxAIaDUu3ao3XGvzD/ZnadCqmxpr7TPP3XsdYCD7vWiOvI1UURclRpsZGzHi9JncDI1l73A+A0lalGVFnBO733TnqfzT9Rs0/gagguLwpb4PVomCPJBH+FC5t1lxXtSqeatX90Pv84vULXR270rRMUwDO+gWx2+sRI9o6UaaopT4iVvKAp6cnQgj10uNLyTutqznQoUZJlh26RWC45jmoITWH4FTUia9Of0V0fJoBjiu1gDL14OQySEzM+4BTKNgJynM1JMRA849TLZZS8vWZrzExMmFCowkAJCZK5uy+Rilbcz5s46SPaJU8cubMGb2NQaleyaPPKHloSjcXIuMS+PHgTQBMjU2Z2nQqD8IfsPpKmkGMhdCcRQXegpv/6CHa/xTcBBUXpbmOWr0LlKiWapX7fXc8Hngwqt4oSlmVAuDvyw+5dD+YiZ2d1Xh7BVh8fHzyvE6KUlhULWnNoKYV2XjmHjeTup03KdOErpW7ssZrDfdC76XewLUX2JaDk0vzPtgUCm6CurwZIgM03wRSiE2IZaHnQqraVeVtl7cBiI5LYOE+H1zL2vKm6lZeoB09elTrLLuKUtB92qEaRcyM+WrPteRlkxpNwtTYlG/Pfpu6sLEpNB0Bfh7w6FIeR/qfgpugTq+E0nXAsVWqxX9c/4MH4Q+Y1HgSpkaajhGrj93hQXAUU7u7qG7lOWTChAnUrFmT0aNHp1sXFRVF27ZtsxzVvE2bNsTHa3mgMBuuXbtGjRr/dYAZNmwYJUuWTJ77Ka379+/j5uaGi4sLrq6uLFq0KNX6tMeS2f5GjBjBkSNHMt3fy2rXrh1+fn56qTuzeHLr56e8uuJJo+K4+zzj6A3NQAkORRx4v/b7uN93x/OxZ+oNGrwLZtaaz1I9KZgJ6ok3PL2qee4pxQ3ZkJgQfr78My3LtaRF2RYAPAuLYbn7LTrVLEWLKiX0FXGB4uvry/Hjx/H29mbJkiXp1q9Zs4Y+ffpgbGyc4T7MzMzo0KEDmzdvztHY0t6gHzp0KPv27cuwvImJCd999x3Xrl3j1KlTLFu2DG/v/2YgTXssme3v9OnTVK1aNdP9ZYc+604rt35+Sva828KRCvaWzNt9jYREzb3AwS6DKVWkFN+d/Y5EmaJThKUduLwB13ZCvJZBZvNAwUxQV7eBMIKaqca0ZcWlFUTERTCh4YTkZYv/vUl0fKLqVp5DfHx8aNu2LXfv3qV+/fpERESkK7NhwwZ69vzvZ+Pm5saBAwcAmDZtGmPGjAGgV69ebNiwIcdiu3r1Kq6urqmWtWnTBnv7jJ93K1OmDA0aaOYNs7GxwcXFhQcP/nvAMe2xZLS/a9euUb16dcqVK5fp/rIjr+tu37499erVo169elhYWLB169ZU63P656dkn7mJMZO71MDnSRh/nvcHwMLEgk8bfMrVwKvsubMn9Qa13oSYULh1UA/RQsHrDSAlXPkTKrcB65LJi++F3mOTzyZ6V+1NtWKaThO3n4Wz8cw93m5SkSoFbRr3vZ/DY6+c3Wfp2tB1fqZFnJ2deffdd3F0dOT9999Ptz42NhZfX18cHR2Tl82aNYsZM2bw9OlTLly4wM6dmoFKatWqhaenZ7p9tG7dmrCwsHTLv/32Wzp27Jj8fuPGjZQsWTJ52fHjxxk+fLhOh6qNn58fFy5coGnTphkeS0b27t1Lly5dMt1fbsmtug8dOgTATz/9hLu7O3369GHZsv8mxsvo56foV/faZVhV4Q7f/3OD1+uUxdLMmO5O3fnd+3cWn19Mx4odsTBJGtTAqS1Y2ms+U2t0z/NYC94Z1KNLEOSryfwp/Hj+R0yNTPm43n9dzr/Z54OFiRFjOlRLuxclG7y8vKhbt67WdQEBAdjZ2aVa1qZNG6SUfP/992zatCn5cpmxsTFmZmbpkpGHhwcXL15M90qZnGJiYrC1tU31DT4hISHTy4qZCQ8P58033+THH3/E1tY2w2PJyP79+1MlCW37e6Fjx47UqlUr3WvHjh2vFPvL1P2y9f/222/s3buXDRs2pGvbjH5+in4JIfiiaw0eh0az5vgdAIyEERMbTeRRxCM2XEtx1mtsqrkS5bMXYtNfDcltBe8M6sqfYGSiGTo+yaVnlzhw9wCj6o3CoYhmkNpzd4PYd/Ux4ztVx8HGXF/R5p4sznRy04tLaQEBAUybNg2Au3fvEh8fz5YtW4iOTv1goJeXF48ePaJEiRLY2KSedysmJgYLi9RDVOlyBmVubk6PHj1YuHAhCQkJBAYGUqpUqVc6nri4ON58800GDRpEnz59kpdbWlqmOxZtIiMjCQ4OpmzZspnu74WDB3PucsrL1v0y9W/dupUNGzawY8cOTE21j/iv7een6F9Tp+J0dCnFT4dvM6BxBYpbm9OkTBPalW/HL16/8Ga1N7GzsNMUrtUHzq2FG/s1/89DBesMKjERrv4FVTqkmvPpp4s/Ucy8GO/WfBfQPKj71Z7rONiY837ryvqKtkAKCwvD1NSUIkWKUKJECVasWMGkSZOws7Pjzz//pFixYiQkJCR/sD969IhBgwaxY8cOrKys2L9/f/K+AgMDcXBwSPfhp8sZ1AvNmzfnxIkTHDhwgE6dOr308UgpGT58OC4uLowfn3qiy7THkhF3d3fc3Nyy3F9uyK26d+3axfLly9m2bVuGCSijn59iGD7v6kxkbDxLDt1KXvZpg0+JiIvgN+/f/itYqSVYl9J8+c9jBStB+XtCyP1Ul/cuP7vM8YfHGVprKEVMiwCw/+oTzt19zvhO1dVDuTnsypUrqbo537x5k+nTp7Nq1arky0mdO3fm2LFjREZG0qdPH7777jtcXFyYPn06X375ZfK27u7udOvWLVvx9OzZkx07dhAWFpbu7Axg4MCBNG/eHB8fH8qXL8/q1f89Vd+tWzc8PDz4/fffOXToUHKHgD17/ruR/OJYMttfyntAx48fz3R/2ZGXdb/77rv4+/vTsmVL6tWrl6rdXsiJn5+Se6qWtKF/44qsP3UXvwDN5buqxarSqVInNl7fSEhMiKagkTG49oabByA6JG+D1NdQJw0bNpQ5bvckKeeUlDIqJHnRyAMjZes/WsuI2AgppZRx8QnS7Vt32eG7wzIuPiHnY9Ajb29vfYeQyrVr1+TAgQNleHh4quXnz5+XgwcPznL73r17y+vXr2crhoSEBNm0aVO5bt26bO0nI7ocS/369WVsbGyu1N+2bVt5584dvdSdVTw58fMzRIb2d5YdT0KjpMv0vfLjDeeSl/kE+cha62rJJeeX/Ffw3mkpZ9pKeWFjrsQBnJVa8kTBOYNKTADv7VCtM1hovql7PfPi2INjDHEdknz2tPWcP77PIvjsNWdMjAvO4Rua0NBQ2rdvj4mJCRMmTGDkyJFERmqG769fvz5ubm5ZPqjbq1cvnJ2dsxWHkZERrq6uyZe5cpoux3L+/Hm9XebSV9059fNTcldJGwuGt6rMrsuPuPJAc3ZUvVh1OlbsyIZrG/47iyrfGIpW1DzCk5e0Za28eOX4GZTvUU2Gv7ItedGog6Nkyz9ayvBYzTf4qNh42WTeAdl72TGZmJiYs/UbgIL0zS4nBQUF6TuEXLN27Vr5/PlzfYeRzNDiyQ0F7e8sJCpW1pu1Xw7+5VTysuuB12WtdbXksgvL/iu4f5qUs+yljAjM8Rgo8GdQD85p/q3SHoCrAVc56n+Ud2u+i5WpFQDrTvjxJDSGyV1qqCH/C5FixYrpO4RcM3ToUJ27uucFQ4tHyZqthSkfu1XF42YAJ24FAOBs70z7Cu1Z772e0NhQTcGqHSAxPuefr8xEwUlQQbfBygEsigLw8+WfsTWzZWCNgQCERMax3P0Wbs4ONHUqntmeFEVRCpXBzSpRtqgFC/b7JE+HMrLuSMLiwth4baOmkH0Vzb9Bt/MsrgKUoO4kN6B/mD+H7x9mYI2BWJtpRohYcfQ2YTHxTHpNDWmkKIqSkoWpMWM7VefS/WD2X30MgEtxF1qVa8UWny3EJcZppt8wNtcMhJBHCk6CCrwN9pqJBrfc2IKRMOKt6m8B8CQ0mrXH79Czbllqlk3/5LyiKEph92aD8lQrac3C/T7EJ2gGjR1YYyDPop5x6N4hMDIC+8oQqBLUy4mNhLCHYO9ETEIMf938i/YV2ydPRrjk0E3iEyTjO6keRYqiKNoYGwkmvuaM77MItl3QDCLcsmxLylmXY9P1TZpC9k7qDOqlPffT/Fvcif1++wmOCaa/c38A7gdFsunMffo3rkDF4kX0F6OiKIqB61yzFHXLF2XRwZvExCdgbGRMP+d+nH1yllvPb2kS1PM7mlF78kDBSFAvbtrZO7Hp+iYqF61Mk9JNAPjx4E2MjQSj26sBYRVFUTIjhOYs6kFwFJvO3Aegd9XemBmZsclnkyZBxUdrrljlgQKSoDSnnFeJxSvAi/7O/RFCcOtpGH9d8GdI80qULqoGrFQURclKq6olaFrZnqXut4iKTaCYRTG6VO7C37f/JqJoOU2hPLrMV3ASVJESbLqzC0sTS96o8gYA3x+4gaWpMR+1q6rnABVFUfIHIQSTXnPmWVgMv570A6C/c38i4yP5OyIpMakE9RICbxNi78jeO3vp4dQDGzMbrjwIYY/XY4a3qoy9lZm+I1QUJRcJIXL0Vdg1crSnnbMDPx2+TWh0HLVL1KZm8ZpsvvsP0thM02s6DxSMBBV0h/3W1sQkxNDPuR8A3/3jQ1FLU95v46Tn4BRFyU23b9/m4sWLOToUmwITOzsTEhXHLx53EELQ37k/t0Ju4128kjqD0llcFIT642mcQEnLkjgXc+b8vee4+zxjRFsnbC3UXDSKUpBdvHgxwxmclVdXq1xRutYqzZpjdwiOjKV1udYAnLWxUwlKZ8/9kMC5uEAalm6IEIIfDtyguJUZ7zZ31Hd0ShZ8fX0ZPnw4ffv21XcoSjZs376dDz74gJ49e/LPP/8kL//111/x8PBg5MiR9O3bl59++inH61ZnPLlnbMfqRMTGs8rDF4ciDlSyrcRZU6EZuScPuprn/wQV5Ms9ExOexYXTqFQjztwJwuNmACPbVsHKXE1GqA8///wzZcqUSZ4U75133smwrJOTk9bJ7qKiomjbtm2WU3K0adOG+Pj4V4ozOjqaJk2aULduXVxdXZk5c2aG9Q8bNoySJUummowxpREjRnDkyBHc3NxwcXHB1dWVRYsWvVJcAO3atcPPzy/5fWb153TdusSVtu179erFqlWrWLduHZs3b04uf+7cOVq1asWKFSvYsmULZ8+eTRf78ePHU7X3/fv3dT6W58+fU7z4f2NrZvVzymrf2f2dKmicS9vQo05Z1h73IygiloalGnI+7jmJ8VEQ/jjX68//CSrwNucszAFoVKoRPxy4gYONOYObVdJzYIXX5cuXmTt3bvJU7L///jteXl706NEj1evp06cZ7mPNmjX06dMHY2PjDMuYmZnRoUOHVB+IL8Pc3JxDhw5x6dIlLl68yL59+zh16pTW+ocOHcq+ffsy3Nfp06epWrUq3333HdeuXePUqVMsW7YMb2/vV4otrczqz+26tcmo7efOncvHH38MQFxcHCYmJggh2LlzJ61ataJDhw7pYm/WrFmq9jYxMdH5WDw8PGjZsmXy+6x+TlntO7u/UwXRpx2qER2XwM9Hb9OoVCNCE2O4aWqaJx0l8n+CCvLlrLUt9hb2PAqw4aRvIKPaVcHSLOMPNiV3eXl5Ua9evVTLateuza5du1K9SpYsmeE+NmzYQM+ePZPfDxgwgP79+9O0aVMqVarE7t27Ac039w0bNrxSnEIIrK01gwnHxcURFxeX3IMrbf1t2rTB3t5e636uXbtG9erVKVeuHA0aNADAxsYGFxcXHjx48EqxpZVR/blVd0btnVLKtpdSMnnyZLp27Zocx9GjR2ndWnPf4o033uDEiROpflYvYjc2Nk7V3mXKlMnwWMLDw5kxY0byPuLi4jAz+6+XbmY/p6z2re24FKha0pqe9crx24m7OFrXBuCsZd4MGpv/r4EF+XLOwpwGJRvw48GblLI1Z2CTivqOSu8WnFnA9aDrObrPGvY1mNxkcpblrl69ynvvvYeRkRElSpTg4MGDGZYNDAxk6tSpXLhwga+//popU6YQGxuLr68vjo6OyeUuXbpEr1692Lx5M8eOHWP8+PF0796dWrVq4enpmW6/rVu3JiwsLN3yb7/9lo4dOya/T0hIoGHDhty6dYuPP/6Ypk2baq0/M3v37qVLly6plvn5+XHhwgWaNm2q0z5eVW7VnVF7p5Sy7ZcsWcLBgwcJCQnh1q1bjBw5kgMHDjBt2jQOHz7Mtm3biImJoVu3buliz6y90x6LpaUllpaW3Lhxg8qVK6dKTi8ro3bK6HeqMBvToRo7Lz1ku2ckZa3KcC4ymkEqQWXtYbAvD+0ELY1r8Jffc+b0dMXCVJ096cv9+/cpXbo0ly9f1ql88eLFWbFiRaplAQEBqSa9i4qKIiAgIPkeUc2aNXn+/DkAxsbGmJmZERYWho2NTfI2Hh4eOtVvbGzMxYsXCQ4Opnfv3ly5cgV7e/uXmnRv//79rF27Nvl9eHg4b775Jj/++CO2tqlHz+/YsSOPH6e/dj9v3rxUZ2z6rDuz9k4pZduPGTOGMWPGpFofHh6OtbU17dq1o127dhnGnvbnndmxGBsb06tXL3bs2EGjRo1SXd57GZm1U0a/U4VZ5RJW9K5fjvWn7tK9fQPOhT1CBt4it58Yy98JKi6ac3FBQHHOeNtRtqgF/RpX0HdUBkGXM53ccPnyZVxdXbO1D0tLS6Kjo5PfX7lyhWrVqmFhoRmu6vz586m6FcfExCSve0HXM6gX7OzsaNeuHfv27WP48OGp6s9MZGQkwcHBlC1bFtBccnrzzTcZNGgQffr0SVc+s7PJl5VbdWfV3ilpa/sXli5dqlPsz58/T9femR2Li4sL33zzDU5OTplezstIVu0EmR9XYTWmfTX+uvCAoKByBBnBneDb5PZTpvk7QQXf5ayFOVbCnCt+RZjTqyrmJursSZ+8vLzSJaiAgACmTZsGwN27d4mPj+fAgQMZ7qNYsWIkJCQQHR2NhYUFly5d4t69e0RHR5OQkMDMmTNZuHAhoLlE6ODggKlp6ufddDmDevbsGaamptjZ2REVFcXBgweZPHlyuvoz4+7ujpubG6C5DzN8+HBcXFwYP358lvVnV27VnVl7p5RR279s7GnbW5djcXBwICYm5qXr1WXf2Tmugqxi8SL0qV+OnVcDMXOEs9FPcJIScnHkjfzdSSKpB59VbDnKFC1Cv0bl9R1Roefl5UXNmjVTLStRogQrVqxg0qRJ2NnZ8eeff2a5n86dO3Ps2DFA84E5aNAg2rVrR+PGjfnoo4+SL+24u7unuq/xMh49eoSbmxt16tShcePGdOrUiR49eqSrH2DgwIE0b94cHx8fypcvn9w1PuU9oOPHj/P7779z6NCh5C72e/bseaXY0tJWf27VnVl7p5Sdtk977yxle+tyLG+99RbNmjVLt9+Mfk7dunXj4cOHOu07O8dV0H3SvirxMfbYSAvN81BhudzVPCeHB3mZV8OGDWV2PT0yX9ZaV0vWXjBW/nbiTrb3l995e3vrO4QM3bhxQw4cOFCGhYXpVP78+fNy8ODBUkopW7duLa9fv661XO/evTNclx0p689M/fr1ZWxsbI7X37ZtW3nnzh291J1Ze6eMKzttnzZ2Xds7L2R1XIb8d5YXJmy5KNss6S7b/1JDJvp65Mg+gbNSS57Q6QxKCNFFCOEjhLglhPhcy3ohhFictP6yEKJBjmdSLTyfXgTAXNRR954M2PXr15k5cyarVq1K7tadlfr16+Pm5kZCQgK3b9+mWrX083nFxsbSq1cvnJ1zfqbklPVn5vz583q7FJRbdWfU3illt+3Txq5re+e23PydKig+catKTERVnpqY4P/4XK7WlWWCEkIYA8uArkBNYKAQomaaYl2BakmvD4GcH89Ei5Mhd7BIhI9btlH3ngxUaGgo7du3x8TEhAkTJjBy5EgiIyN12nbYsGEYGxvz4MEDjIzS/6qamZkxZMiQnA45Xf36MHTo0JfqSZiTMmpv+C+u3Gh7fbb3C7n9O1UQOJawon7F9gAcf3g2i9LZo0sniSbALSmlL4AQYhPQE0j5aHdP4LekU7VTQgg7IUQZKeWjHI84iZSSC4lhVIorwoDGjrlVjZJNtra2PHyYN7NvFiRDhw7VdwhaGWpcSt76rGN7Lv8t8Qi6xYBcrEeXS3zlgPsp3vsnLXvZMgghPhRCnBVCnH327NnLxppKcGQk1tKUelbO6rknRVGUPFTZwZo68cVAmudqPbqcQWnrQ5h2+GBdyiClXAmsBGjUqFG2hiAuZmXFpg8vkpioRjJWFEXJa8veP4qRUe4+qqvLGZQ/kLIHQnkg7TUbXcrkitxuIEVRFCW9vPjs1SVBeQLVhBCVhRBmwABgZ5oyO4EhSb35mgEhuXn/ScmYVHPjKEquUX9feSvLS3xSynghxCfAfsAYWCOlvCqEGJm0fgWwB+gG3AIigfdyL2QlIxYWFgQGBlK8ePHkUbkVRckZUkoCAwPVEEh5SOjrG0GjRo1k2snLlOyJi4vD399f53HkFEV5ORYWFpQvX14Ng5TDhBDnpJSN0i7P32PxKamYmppSuXJlfYehKIqSI/L3WHyKoihKgaUSlKIoimKQVIJSFEVRDJLeOkkIIZ4Bd3NgVyWAgBzYT0Gj2iVjqm20U+2SMdU22uVUu1SSUjqkXai3BJVThBBntfX+KOxUu2RMtY12ql0yptpGu9xuF3WJT1EURTFIKkEpiqIoBqkgJKiV+g7AQKl2yZhqG+1Uu2RMtY12udou+f4elKIoilIwFYQzKEVRFKUAUglKURRFMUj5IkEJIboIIXyEELeEEJ9rWS+EEIuT1l8WQjTQR5z6oEPbDEpqk8tCiBNCiLr6iDOvZdUuKco1FkIkCCH65mV8+qRL2wgh2gkhLgohrgohjuR1jPqgw99SUSHE30KIS0ntUihmbRBCrBFCPBVCXMlgfe59/kopDfqFZoqP24ATYAZcAmqmKdMN2ItmZt9mwGl9x21AbdMCKJb0/66FoW10aZcU5Q6hmS6mr77jNpS2AewAb6Bi0vuS+o7bQNrlC2BB0v8dgCDATN+x50HbtAEaAFcyWJ9rn7/54QyqCXBLSukrpYwFNgE905TpCfwmNU4BdkKIMnkdqB5k2TZSyhNSyudJb0+hme24oNPldwZgNPAn8DQvg9MzXdrmbWCblPIegJSyMLSPLu0iARuhmWzNGk2Cis/bMPOelPIommPNSK59/uaHBFUOuJ/ivX/SspctUxC97HEPR/NNp6DLsl2EEOWA3sCKPIzLEOjyO1MdKCaEOCyEOCeEGJJn0emPLu2yFHABHgJewKdSysS8Cc+g5drnb36YD0rb1LBp+8brUqYg0vm4hRBuaBJUq1yNyDDo0i4/ApOllAmFbPZhXdrGBGgIdAAsgZNCiFNSyhu5HZwe6dIurwEXgfZAFeCAEMJDShmay7EZulz7/M0PCcofqJDifXk032BetkxBpNNxCyHqAL8AXaWUgXkUmz7p0i6NgE1JyakE0E0IES+l3J4nEeqPrn9PAVLKCCBCCHEUqAsU5ASlS7u8B8yXmhsvt4QQd4AawJm8CdFg5drnb364xOcJVBNCVBZCmAEDgJ1pyuwEhiT1JmkGhEgpH+V1oHqQZdsIISoC24B3Cvg34JSybBcpZWUppaOU0hH4HzCqECQn0O3vaQfQWghhIoQoAjQFruVxnHlNl3a5h+asEiFEKcAZ8M3TKA1Trn3+GvwZlJQyXgjxCbAfTU+bNVLKq0KIkUnrV6DphdUNuAVEovmmU+Dp2DYzgOLA8qSzhXhZwEdl1rFdCiVd2kZKeU0IsQ+4DCQCv0gptXYxLih0/J2ZA6wTQnihuaw1WUpZ4KfgEEL8AbQDSggh/IGZgCnk/uevGupIURRFMUj54RKfoiiKUgipBKUoiqIYJJWgFEVRFIOkEpSiKIpikFSCUhRFUQySSlCKoiiKQVIJSlEURTFIKkEpigFJehpf/V0qCipBKUqOEUK4CyE6Jf1/rhBisZYy44UQV5JeY5OWOQohrgkhlgPnST2uGUKITUKIzUKI00KIu0KI7nlwOIqidwY/1JGi5CMzgdlCiJJAfeCNlCuFEA3RDAPTFM1QOaeTZqt9jmZct/eklKO07LcusF1K2V8I0Qr4Htide4ehKIZBnUEpSg5JmthNAOOBAVLKhDRFWgF/SSkjpJThaAbxbZ207m7SZG+pCCEs0Yy2PitpkTdQLDfiVxRDoxKUouQQIURtoAwQI6UM01Ykk80jMlheC7gppYxOet8AzXTkilLgqQSlKDkgaYrrDWimv44QQrympdhRoJcQoogQwgrNjL4eWey6LlBRCGGRtM0s4IccDF1RDJZKUIqSTUlzJm0DJkgpr6GZluHLtOWklOeBdWgmuDuNZhqLC1nsvi6axHcYzZxFP0kpj+dU7IpiyNR0G4piwJJms/1ASumj71gUJa+pBKUoBkwI8QCoIKVM1HcsipLXVIJSFEVRDJK6B6UoiqIYJJWgFEVRFIOkEpSiKIpikFSCUhRFUQySSlCKoiiKQVIJSlEURTFIKkEpiqIoBun/dEtTsMPEvoIAAAAASUVORK5CYII=\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"# Test for the function I found.\n",
|
|
"\n",
|
|
"x = np.linspace(0, 1, 100)\n",
|
|
"g = lambda x: np.sin(np.pi*x*(1 - x))\n",
|
|
"f_Z = lambda z: np.sqrt(1.2*(1/2-np.abs(1/2 - z)))\n",
|
|
"F_Z_inv = lambda p: 3/(2*np.sqrt(1.2))*(1/2 - np.abs(p-.5))**(2/3)\n",
|
|
"\n",
|
|
"plt.figure()\n",
|
|
"plt.title(\"My choice of distribution for $Z$\")\n",
|
|
"plt.plot(x, g(x), label=\"$g(x) = \\sin(\\pi x(1-x))$\")\n",
|
|
"plt.plot(x, f_Z(x), label=\"$f_Z(x) = \\sqrt{1.2(1/2-|1/2 - z|)}$\")\n",
|
|
"plt.plot(x, F_Z_inv(x), label=\"$F_Z^{-1}(p) = 3(1/2 - |1/2 - p|)^{2/3}/(2\\sqrt{1.2})$\")\n",
|
|
"plt.xlabel('$x$ or $p$')\n",
|
|
"plt.legend()\n",
|
|
"plt.tight_layout()\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 134,
|
|
"id": "5793420a",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "d546b8c253155efa96a399ec15c24fb8",
|
|
"grade": true,
|
|
"grade_id": "cell-467d7240beab5b29",
|
|
"locked": false,
|
|
"points": 20,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"\n",
|
|
"After trial and error, I found that 1500 trials were sufficient to arrive at an 1-sigma error of around 0.001,\n",
|
|
"with mean 0.484070 and standard deviation sigma 0.001016. Wolfram Alpha tells me it should be approximately I =\n",
|
|
"0.487595, so I am happy.\n",
|
|
" \n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"def sample_nice_Z():\n",
|
|
" '''Sample from the nice distribution Z'''\n",
|
|
" p = rng.random()\n",
|
|
" #return 3/(2*np.sqrt(1.2))*p**(2/3)\n",
|
|
" #return 3/(2*np.sqrt(1.2))*(p%.5)**(2/3)\n",
|
|
" return 3/(2*np.sqrt(1.2))*np.abs(p - .5)**(2/3)\n",
|
|
" #return np.arctan(np.pi*p)/np.pi\n",
|
|
" #return 1/3*np.arccos(1-4*p)\n",
|
|
" \n",
|
|
"def sample_X_prime():\n",
|
|
" '''Sample from X'.'''\n",
|
|
" z = sample_nice_Z()\n",
|
|
" return np.sin(np.pi*z*(1 - z))/np.sqrt(1.2*(1/2-np.abs(1/2 - z)))/2\n",
|
|
" #return np.sin(np.pi*z*(1 - z))*np.cos(np.pi*z)**2\n",
|
|
" #return np.sin(np.pi*z*(1 - z))/(3/2*np.sin(3/2*z)*np.cos(3/2*z))\n",
|
|
" \n",
|
|
"n = 1500\n",
|
|
"sample_mean, sample_stdev = estimate_expectation(sample_X_prime, n)\n",
|
|
"print(\"\"\"\n",
|
|
"After trial and error, I found that {} trials were sufficient to arrive at an 1-sigma error of around 0.001,\n",
|
|
"with mean {:.6f} and standard deviation sigma {:.6f}. Wolfram Alpha tells me it should be approximately I =\n",
|
|
"{}, so I am happy.\n",
|
|
" \"\"\".format(n, sample_mean, sample_stdev, 0.487595))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "e2a5b934",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "8fd70da242e4db50dd7a7db9851e839a",
|
|
"grade": false,
|
|
"grade_id": "cell-59ae6e6f3f476d2a",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"## Direct sampling of Dyck paths\n",
|
|
"\n",
|
|
"**(35 points)**\n",
|
|
"\n",
|
|
"Direct sampling of random variables in high dimensions requires some luck and/or ingenuity. Here is an example of a probability distribution on $\\mathbb{Z}^{2n+1}$ that features prominently in the combinatorial literature and can be sampled directly in an efficient manner. A sequence $\\mathbf{x}\\equiv(x_0,x_1,\\ldots,x_{2n})\\in\\mathbb{Z}^{2n+1}$ is said to be a **Dyck path** if $x_0=x_{2n}=0$, $x_i \\geq 0$ and $|x_{i}-x_{i-1}|=1$ for all $i=1,\\ldots,2n$. Dyck paths are counted by the Catalan numbers $C(n) = \\frac{1}{n+1}\\binom{2n}{n}$. Let $\\mathbf{X}=(X_0,\\ldots,X_n)$ be a **uniform Dyck path**, i.e. a random variable with probability mass function $p_{\\mathbf{X}}(\\mathbf{x}) = 1/C(n)$ for every Dyck path $\\mathbf{x}$. Here is one way to sample $\\mathbf{X}$."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "d11a0889",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "0ea1008043f2711810b9da3c2cf96a02",
|
|
"grade": false,
|
|
"grade_id": "cell-3598f5c8aabad8b9",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def random_dyck_path(n):\n",
|
|
" '''Returns a uniform Dyck path of length 2n as an array [x_0, x_1, ..., x_{2n}] of length 2n.'''\n",
|
|
" # produce a (2n+1)-step random walk from 0 to -1\n",
|
|
" increments = [1]*n +[-1]*(n+1)\n",
|
|
" rng.shuffle(increments)\n",
|
|
" unconstrained_walk = np.cumsum(increments)\n",
|
|
" # determine the first time it reaches its minimum\n",
|
|
" argmin = np.argmin(unconstrained_walk)\n",
|
|
" # cyclically permute the increments to ensure walk stays non-negative until last step\n",
|
|
" rotated_increments = np.roll(increments,-argmin)\n",
|
|
" # turn off the superfluous -1 step\n",
|
|
" rotated_increments[0] = 0\n",
|
|
" # produce dyck path from increments\n",
|
|
" walk = np.cumsum(rotated_increments)\n",
|
|
" return walk\n",
|
|
"\n",
|
|
"\n",
|
|
"plt.plot(random_dyck_path(50))\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "3e6ad6ad",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "e88c23c1f0a6e2d23e0cd8ab4e2a20fd",
|
|
"grade": false,
|
|
"grade_id": "cell-8c2933e84038cd5e",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"__(a)__ Let $H$ be the (maximal) height of $X$, i.e. $H=\\max_i X_i$. Estimate the expected height $\\mathbb{E}[H]$ for $n = 2^5, 2^6, \\ldots, 2^{11}$ (including error bars). Determine the growth $\\mathbb{E}[H] \\approx a\\,n^\\beta$ via an appropriate fit. *Hint*: use the `scipy.optimize.curve_fit` function with the option `sigma = ...` to incorporate the standard errors on $\\mathbb{E}[H]$ in the fit. Note that when you supply the errors appropriately, fitting on linear or logarithmic scale should result in the same answer. **(25 pts)**"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "2b79f35a",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "54439c27bcec7884b6168c5bf92c27fd",
|
|
"grade": true,
|
|
"grade_id": "cell-310664cb8eaaa831",
|
|
"locked": false,
|
|
"points": 10,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Collect height estimates\n",
|
|
"n_values = [2**k for k in range(5,11+1)]\n",
|
|
"# YOUR CODE HERE\n",
|
|
"raise NotImplementedError()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "77d3723b",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "960be7ade957300785f4b69233a8f462",
|
|
"grade": true,
|
|
"grade_id": "cell-0ca47ced5547d67c",
|
|
"locked": false,
|
|
"points": 10,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"from scipy.optimize import curve_fit\n",
|
|
"\n",
|
|
"# Fitting\n",
|
|
"# YOUR CODE HERE\n",
|
|
"raise NotImplementedError()\n",
|
|
"print(\"Fit parameters: a = {}, beta = {}\".format(a_fit,beta_fit))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "edd467e4",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "421cc39b43f6f9ebc46886d7ca0ab1e4",
|
|
"grade": true,
|
|
"grade_id": "cell-7ebce3629edd06a2",
|
|
"locked": false,
|
|
"points": 5,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Plotting\n",
|
|
"# YOUR CODE HERE\n",
|
|
"raise NotImplementedError()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "142aa5cb",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "1df181e550da9f9c1edff06f45d51984",
|
|
"grade": false,
|
|
"grade_id": "cell-2dcc8f8e438aad86",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"__(b)__ Produce a histogram of the height $H / \\sqrt{n}$ for $n = 2^5, 2^6, \\ldots, 2^{11}$ and $3000$ samples each and demonstrate with a plot that it appears to converge in distribution as $n\\to\\infty$. *Hint*: you could call `plt.hist(...,density=True,histtype='step')` for each $n$ to plot the densities on top of each other. **(10 pts)**"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "eed4339c",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "cb48ed2a17571debd5f778f87f34a2ab",
|
|
"grade": true,
|
|
"grade_id": "cell-b2bb42d9970a7271",
|
|
"locked": false,
|
|
"points": 10,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# YOUR CODE HERE\n",
|
|
"raise NotImplementedError()"
|
|
]
|
|
}
|
|
],
|
|
"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.9.12"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 5
|
|
}
|