979 lines
118 KiB
Plaintext
979 lines
118 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "4c431265",
|
|
"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": "026433a4",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"NAME = \"Kees van Kempen\"\n",
|
|
"NAMES_OF_COLLABORATORS = \"\""
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "3b1bff64",
|
|
"metadata": {},
|
|
"source": [
|
|
"---"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "41d26cde",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "de05c5cadee95d63f1acb0ab3f82894f",
|
|
"grade": false,
|
|
"grade_id": "cell-f29a87a28188c3d0",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"__Exercise sheet 2__\n",
|
|
"\n",
|
|
"Code from the lecture:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"id": "cb41d2a1",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "5435cd2800cbe70e733a364b79e86c9b",
|
|
"grade": false,
|
|
"grade_id": "cell-a6520f459483332d",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"import matplotlib.pylab as plt\n",
|
|
"from scipy.integrate import quad\n",
|
|
"\n",
|
|
"rng = np.random.default_rng()\n",
|
|
"%matplotlib inline\n",
|
|
"\n",
|
|
"def inversion_sample(f_inverse):\n",
|
|
" '''Obtain an inversion sample based on the inverse-CDF f_inverse.'''\n",
|
|
" return f_inverse(rng.random())\n",
|
|
"\n",
|
|
"def compare_plot(samples,pdf,xmin,xmax,bins):\n",
|
|
" '''Draw a plot comparing the histogram of the samples to the expectation coming from the pdf.'''\n",
|
|
" xval = np.linspace(xmin,xmax,bins+1)\n",
|
|
" binsize = (xmax-xmin)/bins\n",
|
|
" # Calculate the expected numbers by numerical integration of the pdf over the bins\n",
|
|
" expected = np.array([quad(pdf,xval[i],xval[i+1])[0] for i in range(bins)])/binsize\n",
|
|
" measured = np.histogram(samples,bins,(xmin,xmax))[0]/(len(samples)*binsize)\n",
|
|
" plt.plot(xval,np.append(expected,expected[-1]),\"-k\",drawstyle=\"steps-post\")\n",
|
|
" plt.bar((xval[:-1]+xval[1:])/2,measured,width=binsize)\n",
|
|
" plt.xlim(xmin,xmax)\n",
|
|
" plt.legend([\"expected\",\"histogram\"])\n",
|
|
" plt.show()\n",
|
|
" \n",
|
|
"def gaussian(x):\n",
|
|
" return np.exp(-x*x/2)/np.sqrt(2*np.pi)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "3317e002",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "d2c3d8374cf18fd1a12c91353f28dbcf",
|
|
"grade": false,
|
|
"grade_id": "cell-e6c28b1e3e8371c3",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"## Sampling random variables via the inversion method \n",
|
|
"__(35 Points)__\n",
|
|
"\n",
|
|
"Recall from the lecture that for any real random variable $X$ we can construct an explicit random variable via the inversion method that is identically distributed. This random variable is given by $F_X^{-1}(U)$ where $F_X$ is the CDF of $X$ and $U$ is a uniform random variable on $(0,1)$ and \n",
|
|
"\n",
|
|
"$$\n",
|
|
"F_X^{-1}(p) := \\inf\\{ x\\in\\mathbb{R} : F_X(x) \\geq p\\}.\n",
|
|
"$$\n",
|
|
"\n",
|
|
"This gives a very general way of sampling $X$ in a computer program, as you will find out in this exercise.\n",
|
|
"\n",
|
|
"__(a)__ Let $X$ be an **exponential random variable** with **rate** $\\lambda$, i.e. a continuous random variable with probability density function $f_X(x) = \\lambda e^{-\\lambda x}$ for $x > 0$. Write a function `f_inverse_exponential` that computes $F_X^{-1}(p)$. Illustrate the corresponding sampler with the help of the function `compare_plot` above. __(10 pts)__"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "6f2c475a",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "4292b1a356454d496a93ef6555f0a7ae",
|
|
"grade": true,
|
|
"grade_id": "cell-311fd25e116f5066",
|
|
"locked": false,
|
|
"points": 5,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"Reasoning from the PDF, we can find the CDF and invert that as follows.\n",
|
|
"\n",
|
|
"$$\n",
|
|
"f_X(x) = \\lambda{}e^{-\\lambda{}x}\n",
|
|
"$$\n",
|
|
"$$\n",
|
|
"\\implies F_X(x)\n",
|
|
" = \\int_{-\\infty}^x f_X(t)dt\n",
|
|
" = \\int_0^x \\lambda{}e^{-\\lambda{}t}dt\n",
|
|
" = \\left[ -e^{\\lambda{}t} \\right]_{t = 0}^x\n",
|
|
" = 1 - e^{\\lambda{}x}\n",
|
|
" = \\mathbb{P}(X \\leq x) = p\n",
|
|
"$$\n",
|
|
"for $x \\in [0, \\infty)$, otherwise zero.\n",
|
|
"\n",
|
|
"Now we seek $x$ as a function of $p$.\n",
|
|
"\n",
|
|
"$$\n",
|
|
"1 - e^{\\lambda{}x} = p\n",
|
|
"\\iff -\\lambda{}x = \\ln{(1-p)}\n",
|
|
"\\iff x = \\frac{\\ln{(1-p)}}{-\\lambda} = F^{-1}_X(p)\n",
|
|
"$$\n",
|
|
"which works, as $1 - p \\geq 0$ as $p \\in [0, 1]$, allowing $\\ln{0} = -\\infty$."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"id": "e6b6428c",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "90de5b60de4e43881ab85442cdff704a",
|
|
"grade": false,
|
|
"grade_id": "cell-06ef7d054d38f5c6",
|
|
"locked": false,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAYh0lEQVR4nO3df3BU9f3v8ec7IRiJ8UcReh2iJjgoRhJzbYiAvbACYqQq4+3tIKFVaSG1VYozrVd0KLbKVBy5LWL1i0yL+G2N5A6i5Srfr0oVKCqF0O43QCg2BdQULREEQRMw8L5/JEIICXuSbLLJyesxk5Fzznt335yR137y2bOfY+6OiIiES1KiGxARkfhTuIuIhJDCXUQkhBTuIiIhpHAXEQmhXol64fPPP98zMzMT9fIiIt3Spk2bPnb3frHqEhbumZmZlJWVJerlRUS6JTN7L0idpmVEREJI4S4iEkIKdxGREErYnLuIhMMXX3xBVVUVtbW1iW4lVFJTU8nIyCAlJaVNj1e4i0i7VFVVkZ6eTmZmJmaW6HZCwd3Zu3cvVVVVZGVltek5NC0jIu1SW1tL3759FexxZGb07du3Xb8NKdxFpN0U7PHX3nOqcBcRCSHNuSfYokWLKCkpCVxfVFREcXFxB3YkIkFEo1F2797N+PHjW/W4SCTCvHnzyM/P76DO6incE6ykpIS168vo3T/2hyZH9uwEULiLdAHRaJSysrJWh3tnUbh3Ab37Z/HfiubGrPuoZGYndCPSPf3+979nwYIFHDlyhKuvvprvfve7TJs2jQ0bNnD06FEKCgooLS3l448/Zvbs2fTt25ft27czcuRInnrqKZKSknjttdd48MEHOXz4MJdccgnPPPMMZ511Fhs3bmTGjBl89tlnnHHGGbz++uvMnj2bmpoa1q1bx/3338+NN97I9OnT2bx5M3V1dfzsZz9jwoQJ1NTUMGXKFCoqKrj88supqanplPOhcBeRuLnnnnuIRqNxfc68vDzmz59/2ppt27ZRWlrKW2+9RUpKCj/84Q/Zvn07N998M7NmzaKmpoZvf/vbDBkyhNWrV7NhwwYqKiq4+OKLKSwsZPny5UQiEebMmcOqVatIS0vj0Ucf5Ze//CUzZ85k4sSJlJaWMnToUD799FP69OnDQw89RFlZGb/+9a8BeOCBBxg9ejSLFy9m//79FBQUMHbsWJ5++mn69OlDeXk55eXlXHXVVXE9Py1RuItIt/fHP/6RTZs2MXToUABqamro378/s2fPZujQoaSmprJgwYLj9QUFBQwcOBCASZMmsW7dOlJTU6moqOCaa64B4MiRIwwfPpzt27dzwQUXHH/us88+u9keXnvtNVasWMG8efOA+ktE33//fdauXcuPfvQjAHJzc8nNze2Yk9CEwl1E4ibWCLujuDu33347jzzyyEn7P/roIw4dOsQXX3xBbW0taWlpwKmXGZoZ7s51113H888/f9Kx8vLyQJclujsvvPACl1122SnHEnGpqC6FFJFub8yYMSxbtow9e/YAsG/fPt577z2Ki4t5+OGHmTx5Mvfdd9/x+g0bNrBz506OHTtGaWkpX//61xk2bBhvvfUWlZWVAHz++ee8++67DB48mN27d7Nx40YADh48SF1dHenp6Rw8ePD4c15//fU88cQTuDsAf/3rXwEYOXIkzz33HABbtmyhvLy8408IGrmLSAhkZ2czZ84cxo0bx7Fjx0hJSWHChAn06tWLoqIijh49yogRI3jjjTdISkpi+PDhzJw5k82bNzNy5EhuueUWkpKSWLJkCZMmTeLw4cMAzJkzh0svvZTS0lKmT59OTU0NZ555JqtWreLaa69l7ty55OXlcf/99/PTn/6Ue+65h9zcXNydzMxMXn75ZX7wgx8wZcoUcnNzycvLo6CgoFPOicJdREJh4sSJTJw4sdljycnJ/PnPfwZg9erV9OnTh9LS0lPqRo8efXyE3tjQoUNZv379Kfub1j799NOn1Jx55pksXbo00N8hnjQtIyISQhq5i0iPEolEiEQiiW6jw2nkLiISQgp3EZEQUriLiISQwl1EJIT0gaqIxFXmzFfi+ny75n4jds2uXdx4441s2bLlpP2zZ89m5MiRjB07ttnHvfTSS1x66aVkZ2fHpdeuRCN3EQmthx56qMVgh/pwr6ioiMtr1dXVxeV54kXhLiKhcPToUaZNm8YVV1zBuHHjqKmp4Y477mDZsmUAzJw5k+zsbHJzc/nJT37C22+/zYoVK7j33nvJy8vjH//4B9FolGHDhpGbm8stt9zCJ598AtR/WSk3N5fhw4dz7733MmTIEACWLFnCt771LW666SbGjRvHoUOHGDNmDFdddRU5OTn84Q9/AOp/sxg8eDBTp05lyJAhTJ48mVWrVnHNNdcwaNAgNmzYEPfzETPczWyxme0xsy0tHJ9sZuUNP2+b2ZVx71JEJIa///3v3HXXXWzdupVzzz2XF1544fixffv28eKLL7J161bKy8uZNWsWI0aM4Oabb+axxx4jGo1yySWXcNttt/Hoo49SXl5OTk4OP//5zwGYMmUKCxcu5J133iE5Ofmk133nnXd49tlneeONN0hNTeXFF1/kL3/5C2+++SY//vGPj681U1lZyYwZMygvL+dvf/sbJSUlrFu3jnnz5vGLX/wi7ucjyMh9CVB4muM7gVHungs8DCyKQ18iIq2SlZVFXl4eAF/72tfYtWvX8WNnn302qampTJ06leXLl9OnT59THn/gwAH279/PqFGjALj99ttZu3Yt+/fv5+DBg4wYMQKov9VlY9dddx1f+cpXgPqVIR944AFyc3MZO3Ys//znP/nXv/51vL+cnBySkpK44oorGDNmDGZGTk7OSb3GS8xwd/e1wL7THH/b3T9p2FwPZMSpNxGRwM4444zjf05OTj5pDrxXr15s2LCBb37zm7z00ksUFp5uvHqyL0feLflyGWGA5557jurqajZt2kQ0GuWrX/0qtbW1p/SXlJR0fDspKalD5uvjPef+PeA/WjpoZsVmVmZmZdXV1XF+aRGR5h06dIgDBw4wfvx45s+ff/xuUY2X7T3nnHM477zz+NOf/gTA7373O0aNGsV5551Henr68YXDTrcI2IEDB+jfvz8pKSm8+eabvPfeex37FzuNuF0KaWbXUh/uX2+pxt0X0TBtk5+ff/q3QxHploJcutjZDh48yIQJE6itrcXd+dWvfgXArbfeyrRp01iwYAHLli3j2Wef5c477+Tzzz9n4MCBPPPMMwD89re/Zdq0aaSlpRGJRDjnnHOafZ3Jkydz0003kZ+fT15eHoMHD+60v2NTFutXDgAzywRedvchLRzPBV4EbnD3d4O8cH5+vpeVlbWi1XCKRCKs37E38A2yhw3sy+rVqzu+MZGAtm3bxuWXX57oNjrUoUOHOOusswCYO3cuH374IY8//niHv25z59bMNrl7fqzHtnvkbmYXAcuB7wQNdhGR7uSVV17hkUceoa6ujosvvpglS5YkuqWYYoa7mT0PRIDzzawKeBBIAXD3hcBsoC/wVMN9AuuCvKuIiHQXp7sRSFcVM9zdfVKM41OBqXHrSES6HXdPyE2gwyzIlPnp6BuqItIuqamp7N27t91hJCe4O3v37iU1NbXNz6GFw0SkXTIyMqiqqkKXN8dXamoqGRlt/9qQwl1E2iUlJYWsrKxEtyFNaFpGRCSEFO4iIiGkcBcRCSGFu4hICCncRURCSOEuIhJCCncRkRBSuIuIhJDCXUQkhBTuIiIhpHAXEQkhhbuISAgp3EVEQkjhLiISQgp3EZEQUriLiISQbtbRzUSjUSKRSKDaoqIiiouLO7YhEemSFO7dSFp2hJzazYFqo9EogMJdpIdSuHcj6XmF7KIwUG3tjpkd3I2IdGUx59zNbLGZ7TGzLS0cNzNbYGaVZlZuZlfFv00REWmNIB+oLoHTDhdvAAY1/BQD/9b+tkREpD1ihru7rwX2naZkAvDvXm89cK6ZXRCvBkVEpPXicSnkAOCDRttVDftOYWbFZlZmZmXV1dVxeGkREWlOPMLdmtnnzRW6+yJ3z3f3/H79+sXhpUVEpDnxCPcq4MJG2xnA7jg8r4iItFE8wn0FcFvDVTPDgAPu/mEcnldERNoo5nXuZvY8EAHON7Mq4EEgBcDdFwIrgfFAJfA5MKWjmhURkWBihru7T4px3IG74taRiIi0mxYOExEJIYW7iEgIKdxFREJI4S4iEkIKdxGREFK4i4iEkMJdRCSEFO4iIiGkcBcRCSGFu4hICCncRURCSOEuIhJCCncRkRBSuIuIhJDCXUQkhBTuIiIhpHAXEQkhhbuISAgp3EVEQkjhLiISQgp3EZEQUriLiIRQoHA3s0Iz225mlWY2s5nj55jZ/zOz/zKzrWY2Jf6tiohIUDHD3cySgSeBG4BsYJKZZTcpuwuocPcrgQjwf8ysd5x7FRGRgHoFqCkAKt19B4CZLQUmABWNahxINzMDzgL2AXVx7lVaae36MlIvyglcv2DWdIqLizuwIxHpLEGmZQYAHzTarmrY19ivgcuB3cBmYIa7H2v6RGZWbGZlZlZWXV3dxpYliLTsCL37ZwWuP7JnJyUlJR3YkYh0piAjd2tmnzfZvh6IAqOBS4DXzexP7v7pSQ9yXwQsAsjPz2/6HBJH6XmFpOcVBq7/qOSUj1JEpBsLMnKvAi5stJ1B/Qi9sSnAcq9XCewEBsenRRERaa0g4b4RGGRmWQ0fkt4KrGhS8z4wBsDMvgpcBuyIZ6MiIhJczGkZd68zs7uBV4FkYLG7bzWzOxuOLwQeBpaY2Wbqp3Huc/ePO7BvERE5jSBz7rj7SmBlk30LG/15NzAuvq2JiEhb6RuqIiIhpHAXEQkhhbuISAgp3EVEQkjhLiISQgp3EZEQUriLiISQwl1EJIQU7iIiIaRwFxEJIYW7iEgIKdxFREJI4S4iEkIKdxGREFK4i4iEkMJdRCSEFO4iIiGkcBcRCSGFu4hICCncRURCKNANsqVniEajRCKRQLVFRUUUFxd3bEMi0mYKdwEgLTtCTu3mQLXRaBRA4S7ShQUKdzMrBB4HkoHfuPvcZmoiwHwgBfjY3UfFrUvpcOl5haye+0Sg2qCjexFJnJjhbmbJwJPAdUAVsNHMVrh7RaOac4GngEJ3f9/M+ndQv9KBMme+Eqjuox17GTawbwd3IyLtEeQD1QKg0t13uPsRYCkwoUlNEbDc3d8HcPc98W1TRERaI0i4DwA+aLRd1bCvsUuB88xstZltMrPbmnsiMys2szIzK6uurm5bxyIiElOQcLdm9nmT7V7A14BvANcDPzWzS095kPsid8939/x+/fq1ulkREQkmyAeqVcCFjbYzgN3N1Hzs7p8Bn5nZWuBK4N24dCkiIq0SZOS+ERhkZllm1hu4FVjRpOYPwP8ws15m1ge4GtgW31ZFRCSomCN3d68zs7uBV6m/FHKxu281szsbji90921m9p9AOXCM+sslt3Rk4yIi0rJA17m7+0pgZZN9C5tsPwY8Fr/WRESkrbS2jIhICCncRURCSOEuIhJCCncRkRBSuIuIhJDCXUQkhBTuIiIhpHAXEQkhhbuISAgp3EVEQkjhLiISQrpBtrTJ2vVlpF6UE6h2wazpupm2SCfTyF1aLS07Qu/+WYFqj+zZSUlJSQd3JCJNaeQurZaeV0h6XmGg2o9KZnZwNyLSHI3cRURCSOEuIhJCCncRkRBSuIuIhJDCXUQkhBTuIiIhpHAXEQkhhbuISAgFCnczKzSz7WZWaWYtfivFzIaa2VEz+1/xa1FERForZribWTLwJHADkA1MMrPsFuoeBV6Nd5MiItI6QUbuBUClu+9w9yPAUmBCM3XTgReAPXHsT0RE2iBIuA8APmi0XdWw7zgzGwDcAiyMX2siItJWQcLdmtnnTbbnA/e5+9HTPpFZsZmVmVlZdXV1wBZFRKS1gqwKWQVc2Gg7A9jdpCYfWGpmAOcD482szt1falzk7ouARQD5+flN3yBERCROgoT7RmCQmWUB/wRuBYoaF7j78cW9zWwJ8HLTYJeeKxqNEolEAtUWFRXpxh4icRAz3N29zszupv4qmGRgsbtvNbM7G45rnl1alJYd4bOK1azfsTdm7ZE9OwEU7iJxEOhmHe6+EljZZF+zoe7ud7S/LQkL3dhDJDH0DVURkRBSuIuIhJDCXUQkhBTuIiIhpHAXEQkhhbuISAgp3EVEQkjhLiISQgp3EZEQUriLiIRQoOUHRDrL2vVlpF6UE6h2wazpWodGpAUauUuXkZYdoXf/rNiF1C8yVlJS0sEdiXRfGrlLl6FFxkTiRyN3EZEQUriLiISQwl1EJIQU7iIiIaRwFxEJIYW7iEgIKdxFREJI17lLtxWNRolEIoFqi4qK9G1W6VEU7tItpWVH+KxiNet37I1Ze2TPTgCFu/QoCnfplvRtVpHTCzTnbmaFZrbdzCrN7JR/KWY22czKG37eNrMr49+qiIgEFTPczSwZeBK4AcgGJplZdpOyncAod88FHgYWxbtREREJLsjIvQCodPcd7n4EWApMaFzg7m+7+ycNm+uBjPi2KSIirREk3AcAHzTarmrY15LvAf/R3AEzKzazMjMrq66uDt6liIi0SpBwt2b2ebOFZtdSH+73NXfc3Re5e7675/fr1y94lyIi0ipBrpapAi5stJ0B7G5aZGa5wG+AG9w99vVpIp2oNdfEg66Ll+4vyMh9IzDIzLLMrDdwK7CicYGZXQQsB77j7u/Gv02RtkvLjpCXlxe4PhqN6i5P0u3FHLm7e52Z3Q28CiQDi919q5nd2XB8ITAb6As8ZWYAde6e33FtiwSXnlfILoJdEw9Qu0PXxUv3F+hLTO6+EljZZN/CRn+eCkyNb2siItJWWjhMRCSEFO4iIiGktWVEmqEVJ6W7U7iLNKEVJyUMFO4iTWjFSQkDzbmLiISQwl1EJIQ0LSPSTvrwVboihbtIO+jDV+mqFO4i7aAPX6Wr0py7iEgIaeQu0onWri8j9aKcQLULZk3XFI60mUbuIp0kLTtC7/5ZgWqP7NmpZYelXTRyF+kkmp+XzqRwF+miNIUj7aFpGZEuSFM40l4auYt0Qa2dwtEXqaQphbtIN6cvUklzFO4i3VxrR/mtmcsHzed3Vwp3kR4kLTvSqvrDH2zh+9//fuA5fU35dB0Kd5EepDWjfICD0f8MPOVz+IMtrFmzRm8EXYTCXURa1Jo3A70RdC2Bwt3MCoHHgWTgN+4+t8lxazg+HvgcuMPd/xLnXkWkC+voN4IfzXki0HMPG9g3UF1bdKc3mZjhbmbJwJPAdUAVsNHMVrh7RaOyG4BBDT9XA//W8F8RkVO05Y0gqCBvGF9qzRvBmjVrWvXbRqIFGbkXAJXuvgPAzJYCE4DG4T4B+Hd3d2C9mZ1rZhe4+4ctPen27dsDX5cbZtFoFM6+MNFtiHRZrf2coDV2taL2K6k5gX/b6AqChPsA4ING21WcOipvrmYAcFK4m1kx8OXvNIfXrFmzpVXdhtWBA+e/9+iNHye6jS7ifEDnop7OxQk6FydcFqQoSLhbM/u8DTW4+yJgEYCZlbl7foDXDz2dixN0Lk7QuThB5+IEMysLUhdkbZkqoPG8QQawuw01IiLSSYKE+0ZgkJllmVlv4FZgRZOaFcBtVm8YcOB08+0iItKxYk7LuHudmd0NvEr9pZCL3X2rmd3ZcHwhsJL6yyArqb8UckqA117U5q7DR+fiBJ2LE3QuTtC5OCHQubD6C1xERCRMtJ67iEgIKdxFREIooeFuZt8ys61mdszMeuRlTmZWaGbbzazSzHrsjTPNbLGZ7TGzHv/dBzO70MzeNLNtDf8+ZiS6p0Qxs1Qz22Bm/9VwLn6e6J4SycySzeyvZvZyrNpEj9y3AP8TWJvgPhKi0dIONwDZwCQzy05sVwmzBOiYryF2P3XAj939cmAYcFcP/v/iMDDa3a8E8oDChivyeqoZwLYghQkNd3ff5u7bE9lDgh1f2sHdjwBfLu3Q47j7WmBfovvoCtz9wy8X3nP3g9T/Yx6Q2K4Sw+sdathMafjpkVeBmFkG8A3gN0HqEz1y7+laWrZBBAAzywT+O/DnBLeSMA1TEVFgD/C6u/fUczEf+N/AsSDFHR7uZrbKzLY089MjR6hNBFq2QXomMzsLeAG4x90/TXQ/ieLuR909j/pvvheY2ZAEt9TpzOxGYI+7bwr6mA6/WYe7j+3o1+jGtGyDNMvMUqgP9ufcfXmi++kK3H2/ma2m/rOZnvbB+zXAzWY2HkgFzjaz37v7t1t6gKZlEivI0g7SwzTc/Oa3wDZ3/2Wi+0kkM+tnZuc2/PlMYCzwt4Q2lQDufr+7Z7h7JvU58cbpgh0SfynkLWZWBQwHXjGzVxPZT2dz9zrgy6UdtgH/1923JrarxDCz54F3gMvMrMrMvpfonhLoGuA7wGgzizb8jE90UwlyAfCmmZVTPxh63d1jXgYoWn5ARCSUNC0jIhJCCncRkRBSuIuIhJDCXUQkhBTuIiIhpHAXEQkhhbuISAj9f3qejdV+/TV7AAAAAElFTkSuQmCC\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"def f_inv_exponential(lam,p):\n",
|
|
" return -np.log(1 - p)/lam\n",
|
|
"\n",
|
|
"f_X = lambda x, lam: lam*np.exp(-lam*x) if x >= 0 else 0\n",
|
|
"\n",
|
|
"# Input parameters as list for flexibility in testing.\n",
|
|
"for lam in [1.5]:\n",
|
|
" pdf = lambda x: f_X(x, lam)\n",
|
|
" samples = [inversion_sample(lambda p: f_inv_exponential(lam, p)) for _ in range(100000)]\n",
|
|
" compare_plot(samples, pdf, -1, 4, 30)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"id": "804aedbf",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "bce45fa412ba32138080832767338e9d",
|
|
"grade": true,
|
|
"grade_id": "cell-2022e00546cf1bb0",
|
|
"locked": true,
|
|
"points": 5,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"from nose.tools import assert_almost_equal\n",
|
|
"assert_almost_equal(f_inv_exponential(1.0,0.6),0.916,delta=0.001)\n",
|
|
"assert_almost_equal(f_inv_exponential(0.3,0.2),0.743,delta=0.001)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "d590b09d",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "08fdb1c6ca42806566800f06d7ffb22b",
|
|
"grade": false,
|
|
"grade_id": "cell-f7e0d9b58c948be5",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"__(b)__ Let now $X$ have the **Pareto distribution** of **shape** $\\alpha > 0$ on $(b,\\infty)$, which has probability density function $f_X(x) = \\alpha b^{\\alpha} x^{-\\alpha-1}$ for $x > b$. Write a function `f_inv_pareto` that computes $F_X^{-1}(p)$. Compare a histogram with a plot of $f_X(x)$ to verify your function numerically. __(10 pts)__"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "47c7a42f",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "1d1fc6a16462f0d238005fdb33a99857",
|
|
"grade": true,
|
|
"grade_id": "cell-199713328dcd510d",
|
|
"locked": false,
|
|
"points": 5,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"$$\n",
|
|
"f_X(x) = \\alpha b^{\\alpha} x^{-\\alpha-1}\n",
|
|
"\\\\\n",
|
|
"\\implies F_X(x) = \\int_{-\\infty}^x f_X(t)dt\n",
|
|
" = \\int_{b}^x \\alpha{}b^\\alpha{}t^{-\\alpha-1}dt\n",
|
|
" = \\alpha{}b^\\alpha{} \\int_{b}^x t^{-\\alpha-1}dt\n",
|
|
" = \\alpha{}b^\\alpha{} \\left[ \\frac{-t^{-\\alpha}}{\\alpha} \\right]_{t = b}^x\n",
|
|
" = b^\\alpha (b^{-\\alpha} - x^{-\\alpha}) = 1 - b^\\alpha x^{-\\alpha} = p,\n",
|
|
"$$\n",
|
|
"for $x > b$, otherwise $F_X(x) = 0$.\n",
|
|
"\n",
|
|
"To find $F_X^{-1}(p)$, we write $p$ as function of $x$.\n",
|
|
"\n",
|
|
"$$\n",
|
|
"p = 1 - b^\\alpha x^{-\\alpha}\n",
|
|
"\\iff b^\\alpha x^{-\\alpha} = 1 - p\n",
|
|
"\\iff x^{-\\alpha} - b^{-\\alpha}(1-p)\n",
|
|
"\\iff x = \\frac{b}{(1-p)^{1/\\alpha}}\n",
|
|
"$$\n",
|
|
"\n",
|
|
"Thus, $F_X^{-1}(p) = \\frac{b}{(1-p)^{1/\\alpha}}$ for $p \\in [0, 1]$."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 5,
|
|
"id": "e177f32d",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "eb07f40a935275cf5883204fc817beaa",
|
|
"grade": false,
|
|
"grade_id": "cell-074f6a1fd6375c22",
|
|
"locked": false,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAbjElEQVR4nO3de3BV5b3/8fcXCI1gvJSLOsZK6GA1QEhpSLk4EBVo9KCM0zrcrNYeTLFe6u/X+is6FqeWU3F0WkU9BxmleFqD+Y0KZRSVWgW8IQS7jYBiKWBN0SYFQUAiBr7nj73hbEOSvULWZu9kfV4ze9hrPc/a+e5ngE/W7Vnm7oiISHR1yXQBIiKSWQoCEZGIUxCIiEScgkBEJOIUBCIiEdct0wU0p3fv3t6vX79MlyEi0mGsW7fuX+7e51i2zcog6NevH9XV1ZkuQ0SkwzCzD451Wx0aEhGJOAWBiEjEKQhERCIuK88RiEjn8MUXX1BbW0tDQ0OmS+k0cnNzyc/PJycnJ7TPVBCISNrU1taSl5dHv379MLNMl9PhuTs7duygtraWgoKC0D5Xh4ZEJG0aGhro1auXQiAkZkavXr1C38NSEIhIWikEwpWO8VQQiIhEnIJARCREsViMZcuWtXm7srKyjN1Iq5PFETR//nxumv1A4P5zb7+RioqKNFYk0nnEYjGqq6u55JJLMl1KYNojiKDKykoO1G0N1PdA3VYqKyvTXJFI+vzhD3+gtLSU4uJifvSjH/Hmm29SVFREQ0MD+/btY+DAgaxfv54VK1YwevRoLr/8cgoLC5kxYwaHDh0CYPny5YwYMYKhQ4dyxRVXsHfvXgDWrl3LyJEjGTJkCKWlpezevZtZs2ZRVVVFcXExVVVV7Nu3jx/+8IcMGzaMb37zm/zxj38EYP/+/UyePJmioiImTZrE/v37MzZG2iOIqO59Czh96pyU/T6unHkcqpEouPnmm4nFYqF+ZnFxMffdd1+L7e+++y5VVVW89tpr5OTk8OMf/5hNmzZx2WWXcfvtt7N//36uvPJKBg0axIoVK1izZg0bN27k7LPPpry8nKeffpqysjJmz57Niy++SM+ePbn77rv5zW9+w8yZM5k0aRJVVVUMGzaMTz/9lB49enDnnXdSXV3Ngw8+CMBtt93GhRdeyIIFC9i1axelpaWMHTuWhx9+mB49elBTU0NNTQ1Dhw4NdWzaQkEgIp3Wn//8Z9atW8ewYcOA+G/hffv2ZdasWQwbNozc3Fzmzp17pH9paSn9+/cHYMqUKbz66qvk5uayceNGRo0aBcCBAwcYMWIEmzZt4owzzjjy2SeddFKzNSxfvpylS5dy7733AvFLav/+97+zatUqbrrpJgCKioooKipKzyAEoCAQkeOitd/c08Xdufrqq7nrrru+tP7jjz9m7969fPHFFzQ0NNCzZ0/g6EszzQx3Z9y4cSxatOhLbTU1NYEu5XR3nnrqKb7xjW8c1ZYtl9bqHIGIdFoXXXQRTz75JHV1dQDs3LmTDz74gIqKCn71q18xbdo0fv7znx/pv2bNGrZu3cqhQ4eoqqri/PPPZ/jw4bz22mts3rwZgM8++4z333+fc889l+3bt7N27VoA9uzZQ2NjI3l5eezZs+fIZ37nO9/hgQcewN0B+Mtf/gLA6NGjefzxxwFYv349NTU16R+QFmiPQEQ6rcLCQmbPns348eM5dOgQOTk5TJw4kW7dujF16lQOHjzIyJEjeemll+jSpQsjRoxg5syZvPPOO0dOHHfp0oWFCxcyZcoUPv/8cwBmz57NOeecQ1VVFTfeeCP79+/nhBNO4MUXX+SCCy5gzpw5FBcXc+utt/KLX/yCm2++maKiItydfv368cwzz3DddddxzTXXUFRURHFxMaWlpRkbJwWBiHRqkyZNYtKkSc22de3alTfffBOAFStW0KNHD6qqqo7qd+GFFx75zT/ZsGHDWL169VHrm/Z9+OGHj+pzwgkn8MQTTwT6DumWMgjMbAEwAahz90HNtN8CTEv6vPOAPu6+08y2AXuAg0Cju5eEVbiIiIQjyDmChUB5S43ufo+7F7t7MXArsNLddyZ1uSDRrhAQkaxVVlbGM888k+kyMiJlELj7KmBnqn4JU4BFKXuJiEjWCO2qITPrQXzP4amk1Q4sN7N1ZtbqHAVmVmFm1WZWXV9fH1ZZIiKSQpiXj14KvNbksNAodx8KXAxcb2ajW9rY3ee7e4m7l/Tp0yfEskREpDVhBsFkmhwWcvftiT/rgMVA5q6PEhGRZoVy+aiZnQyMAa5MWtcT6OLuexLvxwN3hvHzRKRj6jfz2VA/b9ucf0vdZ9s2JkyYwPr167+0ftasWYwePZqxY8c2u92SJUs455xzKCwsDKXWbBbk8tFFQBnQ28xqgTuAHAB3n5fodjmw3N33JW16GrA4cQt1N6DS3Z8Pr3QRkWN3552t/166ZMkSJkyYEEoQNDY20q1b9t62FeSqoSnufoa757h7vrs/6u7zkkIAd1/o7pObbLfF3YckXgPd/T/S8QVERFI5ePAg1157LQMHDmT8+PHs37+fH/zgBzz55JMAzJw5k8LCQoqKivjZz37G66+/ztKlS7nlllsoLi7mb3/7G7FYjOHDh1NUVMTll1/OJ598AsRvHisqKmLEiBHccsstDBoUv91q4cKFXHHFFVx66aWMHz+evXv3ctFFFzF06FAGDx58ZDrqbdu2ce655zJ9+nQGDRrEtGnTePHFFxk1ahQDBgxgzZo1aR8fzTUkIp3eX//6V66//no2bNjAKaecwlNP/e/FjTt37mTx4sVs2LCBmpoabr/9dkaOHMlll13GPffcQywW4+tf/zpXXXUVd999NzU1NQwePJhf/vKXAFxzzTXMmzePN954g65du37p577xxhs89thjvPTSS+Tm5rJ48WLeeustXn75ZX76058emX9o8+bN/OQnP6Gmpob33nuPyspKXn31Ve69915+/etfp318FAQi0ukVFBRQXFwMwLe+9S22bdt2pO2kk04iNzeX6dOn8/TTT9OjR4+jtt+9eze7du1izJgxAFx99dWsWrWKXbt2sWfPHkaOHAnA1KlTv7TduHHj+OpXvwrEZyG97bbbKCoqYuzYsfzjH//gn//855H6Bg8eTJcuXRg4cCAXXXQRZsbgwYO/VGu6KAhEpNP7yle+cuR9165daWxsPLLcrVs31qxZw3e/+12WLFlCeXmLEykc5fBv9C05PL01wOOPP059fT3r1q0jFotx2mmn0dDQcFR9Xbp0ObLcpUuXL9WaLgoCEYm0vXv3snv3bi655BLuu+++I09RS55O+uSTT+bUU0/llVdeAeD3v/89Y8aM4dRTTyUvL+/IxHOtTSK3e/du+vbtS05ODi+//DIffPBBer9YG2TvaWwR6XSCXO55vO3Zs4eJEyfS0NCAu/Pb3/4WgMmTJ3Pttdcyd+5cnnzySR577DFmzJjBZ599Rv/+/fnd734HwKOPPsq1115Lz549KSsr4+STT27250ybNo1LL72UkpISiouLOffcc4/bd0zFUu3aZEJJSYlXV1dnuoxOq6ysjNVbdgR+ZvHw/r1YsWJF+guTTufdd9/lvPPOy3QZabV3715OPPFEAObMmcNHH33E/fffn9af2dy4mtm6Y53cU3sEIiLt8Oyzz3LXXXfR2NjI2WefzcKFCzNdUpspCERE2qG1B990FDpZLCJplY2HnzuydIyngkBE0iY3N5cdO3YoDELi7uzYsYPc3NxQP1eHhkQkbfLz86mtrUXPGAlPbm4u+fn5oX6mgkBE0iYnJ4eCgoJMlyEp6NCQiEjEKQhERCJOQSAiEnEKAhGRiFMQiIhEnIJARCTiFAQiIhGXMgjMbIGZ1ZnZ+hbay8xst5nFEq9ZSW3lZrbJzDab2cwwCxcRkXAE2SNYCKR6ZM8r7l6ceN0JYGZdgYeAi4FCYIqZFbanWBERCV/KIHD3VcDOY/jsUmCzu29x9wPAE8DEY/gcERFJo7DOEYwws7fN7DkzG5hYdybwYVKf2sS6ZplZhZlVm1m15iURETl+wgiCt4Cz3X0I8ACwJLHemunb4hSE7j7f3UvcvaRPnz4hlCUiIkG0Owjc/VN335t4vwzIMbPexPcAzkrqmg9sb+/PExGRcLU7CMzsdDOzxPvSxGfuANYCA8yswMy6A5OBpe39eSIiEq6U01Cb2SKgDOhtZrXAHUAOgLvPA74HXGdmjcB+YLLHn0LRaGY3AC8AXYEF7r4hLd9CRESOWcogcPcpKdofBB5soW0ZsOzYShMRkeNBdxaLiEScgkBEJOIUBCIiEacgEBGJOAWBiEjEKQhERCJOQSAiEnEKAhGRiFMQiIhEnIJARCTiFAQiIhGnIBARiTgFgYhIxKWcfVQ6hvnz53PT7AcC9T1Qt5XufQvSXJGIdBTaI+gkKisrOVC3NVDf7n0L6FlYlt6CRKTD0B5BJ9K9bwGnT52T6TJEpIPRHoGISMQpCEREIk5BICIScSmDwMwWmFmdma1voX2amdUkXq+b2ZCktm1m9o6ZxcysOszCRUQkHEH2CBYC5a20bwXGuHsR8CtgfpP2C9y92N1Ljq1EERFJp5RXDbn7KjPr10r760mLq4H8EOoSEZHjJOxzBP8OPJe07MByM1tnZhWtbWhmFWZWbWbV9fX1IZclIiItCe0+AjO7gHgQnJ+0epS7bzezvsCfzOw9d1/V3PbuPp/EYaWSkhIPqy4REWldKHsEZlYEPAJMdPcdh9e7+/bEn3XAYqA0jJ8nIiLhaXcQmNnXgKeB77v7+0nre5pZ3uH3wHig2SuPREQkc1IeGjKzRUAZ0NvMaoE7gBwAd58HzAJ6Af9pZgCNiSuETgMWJ9Z1Ayrd/fk0fAcREWmHIFcNTUnRPh2Y3sz6LcCQo7cQEZFsojuLRUQiTrOPSkqxWIyysrJAfadOnUpFRatXCotIllEQSKt6FpYxuOGdQH1jsRiAgkCkg1EQSKvyistZMSfYk8+C7jWISHZREEhK/WY+G6jfx1t2MLx/rzRXIyJh08liEZGIUxCIiEScgkBEJOIUBCIiEacgEBGJOAWBiEjEKQhERCJOQSAiEnEKAhGRiFMQiIhEnIJARCTiFAQiIhGnIBARiTgFgYhIxKUMAjNbYGZ1Zra+hXYzs7lmttnMasxsaFJbuZltSrTNDLNwEREJR5A9goVAeSvtFwMDEq8K4L8AzKwr8FCivRCYYmaF7SlWRETClzII3H0VsLOVLhOB//a41cApZnYGUApsdvct7n4AeCLRV0REskgY5wjOBD5MWq5NrGtpfbPMrMLMqs2sur6+PoSyREQkiDCCwJpZ562sb5a7z3f3Encv6dOnTwhliYhIEGE8s7gWOCtpOR/YDnRvYb2IiGSRMPYIlgJXJa4eGg7sdvePgLXAADMrMLPuwOREXxERySIp9wjMbBFQBvQ2s1rgDiAHwN3nAcuAS4DNwGfANYm2RjO7AXgB6AoscPcNafgOIiLSDimDwN2npGh34PoW2pYRDwoREclSurNYRCTiFAQiIhGnIBARiTgFgYhIxCkIREQiTkEgIhJxCgIRkYhTEIiIRJyCQEQk4hQEIiIRpyAQEYk4BYGISMQpCEREIk5BICIScQoCEZGIUxCIiEScgkBEJOLCeHi9yBGxWIyysrJAfadOnUpFRUV6CxKRlBQEEpqehWXs27iC1Vt2pOx7oG4rgIJAJAsECgIzKwfuJ/4Q+kfcfU6T9luAaUmfeR7Qx913mtk2YA9wEGh095KQapcsk1dcTl5xeaC+H1fOTHM1IhJUyiAws67AQ8A4oBZYa2ZL3X3j4T7ufg9wT6L/pcD/cfedSR9zgbv/K9TKRUQkFEFOFpcCm919i7sfAJ4AJrbSfwqwKIziREQk/YIEwZnAh0nLtYl1RzGzHkA58FTSageWm9k6M2vxgLCZVZhZtZlV19fXByhLRETCECQIrJl13kLfS4HXmhwWGuXuQ4GLgevNbHRzG7r7fHcvcfeSPn36BChLRETCECQIaoGzkpbzge0t9J1Mk8NC7r498WcdsJj4oSYREckSQYJgLTDAzArMrDvx/+yXNu1kZicDY4A/Jq3raWZ5h98D44H1YRQuIiLhSHnVkLs3mtkNwAvELx9d4O4bzGxGon1eouvlwHJ335e0+WnAYjM7/LMq3f35ML+AiIi0T6D7CNx9GbCsybp5TZYXAgubrNsCDGlXhSIiklaaa0hEJOIUBCIiEacgEBGJOAWBiEjEKQhERCJOQSAiEnEKAhGRiFMQiIhEnIJARCTiFAQiIhGnIBARiTgFgYhIxCkIREQiTkEgIhJxCgIRkYhTEIiIRFygB9OIpMOq1dXkfm1woL5zb7+RioqKNFckEk3aI5CM6FlYRve+BYH6HqjbSmVlZZorEoku7RFIRuQVl5NXXB6o78eVM9NcjUi0BdojMLNyM9tkZpvN7Kh/lWZWZma7zSyWeM0Kuq2IiGRWyj0CM+sKPASMA2qBtWa21N03Nun6irtPOMZtRUQkQ4LsEZQCm919i7sfAJ4AJgb8/PZsKyIix0GQIDgT+DBpuTaxrqkRZva2mT1nZgPbuC1mVmFm1WZWXV9fH6AsEREJQ5AgsGbWeZPlt4Cz3X0I8ACwpA3bxle6z3f3Encv6dOnT4CyREQkDEGCoBY4K2k5H9ie3MHdP3X3vYn3y4AcM+sdZFsREcmsIEGwFhhgZgVm1h2YDCxN7mBmp5uZJd6XJj53R5BtRUQks1JeNeTujWZ2A/AC0BVY4O4bzGxGon0e8D3gOjNrBPYDk93dgWa3TdN3ERGRYxDohrLE4Z5lTdbNS3r/IPBg0G1FRCR7aIoJEZGIUxCIiEScgkBEJOIUBCIiEacgEBGJOE1DLR1CLBajrKwscP+pU6fqQTYiASkIJOv1LCxjcMM7gfvHYjEABYFIQAoCyXp5xeVsI9hDbAAatuixFyJtoXMEIiIRpyAQEYk4BYGISMQpCEREIk5BICIScQoCEZGIUxCIiEScgkBEJOIUBCIiEac7i6VTasvcRJqXSKJOQSCdTlvmJtK8RCIBg8DMyoH7iT+A/hF3n9OkfRrw88TiXuA6d3870bYN2AMcBBrdvSSc0kWa15a5iTQvkUiAIDCzrsBDwDigFlhrZkvdfWNSt63AGHf/xMwuBuYD305qv8Dd/xVi3SIiEpIgJ4tLgc3uvsXdDwBPABOTO7j76+7+SWJxNZAfbpkiIpIuQYLgTODDpOXaxLqW/DvwXNKyA8vNbJ2ZtXgg1swqzKzazKrr6+sDlCUiImEIco7AmlnnzXY0u4B4EJyftHqUu283s77An8zsPXdfddQHus8nfkiJkpKSZj9fRETCF2SPoBY4K2k5H9jetJOZFQGPABPdfcfh9e6+PfFnHbCY+KEmERHJEkH2CNYCA8ysAPgHMBmYmtzBzL4GPA18393fT1rfE+ji7nsS78cDd4ZVvEgYdM+BRF3KIHD3RjO7AXiB+OWjC9x9g5nNSLTPA2YBvYD/NDP438tETwMWJ9Z1Ayrd/fm0fBORY9CzsIx9G1ewesuOlH0P1G0FdM+BdD6B7iNw92XAsibr5iW9nw5Mb2a7LcCQdtYokjZ5xeXkFQe75+DjSt1zIJ2T5hoSEYk4BYGISMRpriGRNtCJZemMFAQiAenEsnRWCgKRgHRiWTornSMQEYk47RGIpInOJ0hHoSAQSQOdT5COREEgkgY6nyAdiYJAJAvoMJJkkoJAJMP0jGXJNAWBSIbpGcuSaQoCkQ6mLYeRQIeSJDUFgUgH0parkQA+/3A9K1eupLKyMlB/hUY0KQhEOpC2XI0EsCf2vC5jlZQUBCKdWFsvY9XVS9GkIBARoG2HnXTIqXNREIgI0La9h7YccmpraLSFAiYcCgIRabN0hQbA8P69AvVbuXJl2gIGohUygYLAzMqB+4k/vP4Rd5/TpN0S7ZcAnwE/cPe3gmwrIp1bW09wbwvY76u5g9sUMG2Rzr2YbJQyCMysK/AQMA6oBdaa2VJ335jU7WJgQOL1beC/gG8H3FZEpM3aGjBt0da9mI4uyB5BKbDZ3bcAmNkTwEQg+T/zicB/u7sDq83sFDM7A+gXYNujbNq0qU03zEhi6oGTzsp0GSKdQjpDJl0+uHvCMW8bJAjOBD5MWq4l/lt/qj5nBtwWADOrAA4fkPt85cqV6wPUlkm9gX9luogv2b27ub8M2Vdn81RnuFRnuDpCnd841g2DBIE1s84D9gmybXyl+3xgPoCZVbt7SYDaMqYj1AiqM2yqM1yqMzxmVn2s2wYJglog+ZhDPrA9YJ/uAbYVEZEMCvLM4rXAADMrMLPuwGRgaZM+S4GrLG44sNvdPwq4rYiIZFDKPQJ3bzSzG4AXiF8CusDdN5jZjET7PGAZ8UtHNxO/fPSa1rYNUNf8Y/kyx1lHqBFUZ9hUZ7hUZ3iOuUaLX+gjIiJRFeTQkIiIdGIKAhGRiMtYEJhZuZltMrPNZnbU8/cSJ57nJtprzGxoltZZZma7zSyWeM3KQI0LzKzOzJq99yKLxjJVnRkfy0QdZ5nZy2b2rpltMLOfNNMn42MasM6MjqmZ5ZrZGjN7O1HjL5vpkw1jGaTOrPj7mailq5n9xcyeaaat7ePp7sf9RfzE8d+A/sQvMX0bKGzS5xLgOeL3IgwH3szSOsuAZzIxjkk1jAaGAutbaM/4WAasM+NjmajjDGBo4n0e8H6W/v0MUmdGxzQxPicm3ucAbwLDs3Asg9SZFX8/E7X8X6CyuXqOZTwztUdwZNoKdz8AHJ56ItmRaSvcfTVweNqKbKsz49x9FbCzlS7ZMJZB6swK7v6RJyZNdPc9wLvE75JPlvExDVhnRiXGZ29iMSfxanqFSjaMZZA6s4KZ5QP/BjzSQpc2j2emgqClKSna2ifdgtYwIrFL+ZyZDTw+pbVJNoxlUFk1lmbWD/gm8d8Qk2XVmLZSJ2R4TBOHMWJAHfAnd8/KsQxQJ2TH38/7gP8HHGqhvc3jmakgaM+0FcdTkBreAs529yHAA8CSdBd1DLJhLIPIqrE0sxOBp4Cb3f3Tps3NbJKRMU1RZ8bH1N0Punsx8ZkFSs1sUJMuWTGWAerM+Fia2QSgzt3XtdatmXWtjmemgqA901YcTylrcPdPD+9SuvsyIMfMeh+/EgPJhrFMKZvG0sxyiP/n+ri7P91Ml6wY01R1ZtOYuvsuYAXQdFrPrBjLw1qqM0vGchRwmZltI36o+kIz+0OTPm0ez0wFQXumrciqOs3sdDOzxPtS4mOabZOYZ8NYppQtY5mo4VHgXXf/TQvdMj6mQerM9JiaWR8zOyXx/gRgLPBek27ZMJYp68z0WAK4+63unu/u/Yj/f/SSu1/ZpFubxzMjj6r0dkxbkYV1fg+4zswagf3AZE+cuj9ezGwR8SsaeptZLXAH8ZNdWTOWAevM+FgmjAK+D7yTOGYMcBvwtaRas2FMg9SZ6TE9A3jM4g+p6gL8f3d/Jtv+rQesM9Nj2aL2jqemmBARiTjdWSwiEnEKAhGRiFMQiIhEnIJARCTiFAQiIhGnIBARiTgFgYhIxP0PM2orgbpvoW8AAAAASUVORK5CYII=\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"### Solution\n",
|
|
"def f_inv_pareto(alpha,b,p):\n",
|
|
" return b/(1-p)**(1/alpha)\n",
|
|
"\n",
|
|
"# plotting\n",
|
|
"f_X = lambda alpha, b, x: alpha*b**alpha*x**(-alpha-1) if x >= b else 0\n",
|
|
"\n",
|
|
"# Input parameters as list for flexibility in testing.\n",
|
|
"for params in [(3., 1.)]:\n",
|
|
" alpha, b = params\n",
|
|
" pdf = lambda x: f_X(alpha, b, x)\n",
|
|
" samples = [inversion_sample(lambda p: f_inv_pareto(alpha, b, p)) for _ in range(100000)]\n",
|
|
" compare_plot(samples, pdf, b-1, 4, 30)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 6,
|
|
"id": "c0e1426f",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "62920089752d067b0945eb1d6d98135f",
|
|
"grade": true,
|
|
"grade_id": "cell-726b321246679d28",
|
|
"locked": true,
|
|
"points": 5,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"from nose.tools import assert_almost_equal\n",
|
|
"assert_almost_equal(f_inv_pareto(1.0,1.5,0.6),3.75,delta=0.0001)\n",
|
|
"assert_almost_equal(f_inv_pareto(2.0,2.25,0.3),2.689,delta=0.001)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "66d91446",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "0f3c9abbe9fe756c5cf4bdd6a8a37ac2",
|
|
"grade": false,
|
|
"grade_id": "cell-50306550727804ca",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"__(c)__ Let $X$ be a discrete random variable taking values in $\\{1,2,\\ldots,n\\}$. Write a Python function `f_inv_discrete` that takes the probability mass function $p_X$ as a list `prob_list` given by $[p_X(1),\\ldots,p_X(n)]$ and returns a random sample with the distribution of $X$ using the inversion method. Verify the working of your function numerically on an example. __(15 pts)__"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 7,
|
|
"id": "210f1302",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "93d51c9c889dd5ba3490e0ee298d4240",
|
|
"grade": false,
|
|
"grade_id": "cell-694eb1261c2dc217",
|
|
"locked": false,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAa30lEQVR4nO3df3BV5b3v8feHCKYiqEVRh2BBD1YjxFwa48/BWJWiVXOcng6ordYeydDWX3Nbb9GxtFVvK1OmtbT20EyLelopnFGxjNL6o+rgz0NCmxMFxaaINUULYkWoQQx+7x/ZcDdhh6yEhLVhfV4zGdZaz/Ps/d0PTD6sZ++9liICMzPLrgFpF2BmZulyEJiZZZyDwMws4xwEZmYZ5yAwM8u4/dIuoJBDDz00Ro0alXYZZmZ7jWXLlr0dEYf1ZmxRBsGoUaNobGxMuwwzs72GpNd7O9ZLQ2ZmGecgMDPLOAeBmVnGFeV7BIV8+OGHtLa2snnz5rRL2WeUlpZSVlbGwIED0y7FzFK01wRBa2srQ4YMYdSoUUhKu5y9XkSwfv16WltbGT16dNrlmFmK9pqloc2bNzNs2DCHQB+RxLBhw3yGZWbJgkDSJEkrJbVImr6LfidJ2irp33o6NmEduzPcOvF8mhkkCAJJJcCdwHlAOXCJpPIu+s0EHunpWDMzS0+S9wiqgZaIWAUgaT5QC6zo1O8a4H7gpF6MtQKamppYs2YN559/fo/G1dTUMGvWLKqqqvqpMrPdV19fz7x589Iuw0gWBCOAN/L2W4GT8ztIGgFcDHyaHYOg27F5j1EH1AEcddRRCcra9zU1NdHY2NjjIDDbG8ybN48lLzQyaLg/rJC2JEFQaCG5823N7gC+GRFbO607JxnbcTCiHqgHqKqqKsrbpv36179m9uzZbNmyhZNPPpkvf/nLTJ06laVLl7J161aqq6tZsGABb7/9NjNmzGDYsGGsXLmSCRMm8LOf/YwBAwbw6KOP8u1vf5sPPviAY445hrvuuosDDzyQhoYGrrvuOv75z3+y//7789hjjzFjxgza2tp45plnuPHGG7ngggu45pprePHFF2lvb+c73/kOtbW1tLW1ceWVV7JixQqOP/542tra0p4qs0QGDR/NEZfennYZ+4TXZ17Q67FJgqAVGJm3Xwas6dSnCpifC4FDgfMltScc22PXX389TU1Nu/swO6isrOSOO+7osv3ll19mwYIFPPvsswwcOJCvfvWrrFy5kosuuoibb76ZtrY2vvCFLzB27Fieeuopli5dyooVK/jEJz7BpEmTeOCBB6ipqeG2227j8ccfZ/DgwcycOZMf/vCHTJ8+ncmTJ7NgwQJOOukk3nvvPQ444ABuueUWGhsb+elPfwrATTfdxKc//Wnmzp3Lu+++S3V1Neeccw4///nPOeCAA2hubqa5uZnx48f36dyY2b4tSRA0AGMkjQb+BkwBLs3vEBHbz+0k3Q08FBEPStqvu7F7iz/84Q8sW7aMk07qWPlqa2tj+PDhzJgxg5NOOonS0lJmz569vX91dTVHH300AJdccgnPPPMMpaWlrFixgtNPPx2ALVu2cOqpp7Jy5UqOPPLI7Y89dOjQgjU8+uijLFq0iFmzZgEdH6n961//ypIlS7j22msBqKiooKKion8mwcz2Sd0GQUS0S7qajk8DlQBzI2K5pGm59jk9Hbu7Re/qf+79JSK44oor+P73v7/D8bfeeotNmzbx4YcfsnnzZgYPHgzs/NFMSUQE5557Lr/5zW92aGtubk70Uc6I4P777+eTn/zkTm3+KKiZ9Vai7xFExOKIODYijomI/5s7NqdQCETElyLivl2N3RudffbZ3HfffaxduxaAd955h9dff526ujpuvfVWLrvsMr75zW9u77906VJee+01PvroIxYsWMAZZ5zBKaecwrPPPktLSwsA77//Pq+++irHHXcca9asoaGhAYCNGzfS3t7OkCFD2Lhx4/bH/MxnPsNPfvITIjreQvnTn/4EwIQJE7j33nsBeOmll2hubu7/CTGzfcZec4mJtJWXl3PbbbcxceJEPvroIwYOHEhtbS377bcfl156KVu3buW0007jiSeeYMCAAZx66qlMnz6dF198kQkTJnDxxRczYMAA7r77bi655BI++OADAG677TaOPfZYFixYwDXXXENbWxsf+9jHePzxxznrrLO4/fbbqays5MYbb+Rb3/oW119/PRUVFUQEo0aN4qGHHuIrX/kKV155JRUVFVRWVlJdXZ3ybJnZ3kTb/ndZTKqqqqLzjWlefvlljj/++JQq6pmnnnqKWbNm8dBDD6VdSrf2pnm1fUtNTQ0vrFrvTw31kddnXrAsInr15aG95lpDZmbWP7w01A9qamqoqalJuwwzs0R8RmBmlnEOAjOzjHMQmJllnIPAzCzj9to3i0dNf7hPH2/17Z/ddfvq1VxwwQW89NJLOxyfMWMGEyZM4Jxzzik47sEHH+TYY4+lvNy3YTCz4uQzgt10yy23dBkC0BEEK1b0ze0X2tvb++RxzMzyOQh6YOvWrUydOpUTTjiBiRMn0tbWxpe+9CXuu6/jihrTp0+nvLyciooKvvGNb/Dcc8+xaNEibrjhBiorK/nLX/5CU1MTp5xyChUVFVx88cX84x//AKChoYGKigpOPfVUbrjhBsaOHQvA3Xffzec//3kuvPBCJk6cyKZNmzj77LMZP34848aN47e//S3QccZy3HHHcdVVVzF27Fguu+wyHn/8cU4//XTGjBnD0qVL05k0Myt6DoIe+POf/8zXvvY1li9fzsEHH8z999+/ve2dd95h4cKFLF++nObmZm6++WZOO+00LrroIn7wgx/Q1NTEMcccw+WXX87MmTNpbm5m3LhxfPe73wXgyiuvZM6cOTz//POUlJTs8LzPP/8899xzD0888QSlpaUsXLiQP/7xjzz55JN8/etf337toZaWFq677jqam5t55ZVXmDdvHs888wyzZs3ie9/73p6bKDPbqzgIemD06NFUVlYC8KlPfYrVq1dvbxs6dCilpaVcddVVPPDAAxxwwAE7jd+wYQPvvvsuZ555JgBXXHEFS5Ys4d1332Xjxo2cdtppAFx66Y5X6j733HP5+Mc/DnRcgfSmm26ioqKCc845h7/97W/8/e9/317fuHHjGDBgACeccAJnn302khg3btwOtZqZ5XMQ9MD++++/fbukpGSHNfv99tuPpUuX8rnPfY4HH3yQSZMmJX7c7q73tO3S1gD33nsv69atY9myZTQ1NXH44YezefPmneobMGDA9v0BAwb4/QUz65KDoI9s2rSJDRs2cP7553PHHXdsv4Na/qWkDzroIA455BCefvppAH71q19x5plncsghhzBkyBBeeOEFAObPn9/l82zYsIHhw4czcOBAnnzySV5//fX+fWFmts/baz8+2t3HPfe0jRs3Ultby+bNm4kIfvSjHwEwZcoUpk6dyuzZs7nvvvu45557mDZtGu+//z5HH300d911FwC//OUvmTp1KoMHD6ampoaDDjqo4PNcdtllXHjhhVRVVVFZWclxxx23x16jme2bfBnqIrFp0yYOPPBAAG6//XbefPNNfvzjH/f78+7r82rFy5eh7lv9fhlqSZMkrZTUIml6gfZaSc2SmiQ1Sjojr221pBe3tfWmyCx4+OGHqaysZOzYsTz99NPcfPPNaZdkZhnR7dKQpBLgTuBcoBVokLQoIvK/JfUHYFFEhKQK4L+A/DWLsyLi7T6se58zefJkJk+enHYZZpZBSc4IqoGWiFgVEVuA+UBtfoeI2BT/f41pMNAv603FuIy1N/N8mhkkC4IRwBt5+625YzuQdLGkV4CHgS/nNQXwqKRlkuq6ehJJdbllpcZ169bt1F5aWsr69ev9y6uPRATr16+ntLQ07VLMLGVJPjWkAsd2+m0cEQuBhZImALcC2y7Ac3pErJE0HHhM0isRsaTA+HqgHjreLO7cXlZWRmtrK4VCwnqntLSUsrKytMsws5QlCYJWYGTefhmwpqvOEbFE0jGSDo2ItyNiTe74WkkL6Vhq2ikIujNw4EBGjx7d02FmZtaNJEtDDcAYSaMlDQKmAIvyO0j6F0nKbY8HBgHrJQ2WNCR3fDAwEdjxOs5mZpaqbs8IIqJd0tXAI0AJMDcilkualmufA3wOuFzSh0AbMDn3CaLD6Vgu2vZc8yLi9/30WszMrBcSfbM4IhYDizsdm5O3PROYWWDcKuDE3azRzMz6ka81ZGaWcQ4CM7OMcxCYmWWcg8DMLOMcBGZmGecgMDPLOAeBmVnGOQjMzDLOQWBmlnEOAjOzjHMQmJllnIPAzCzjHARmZhnnIDAzyzgHgZlZxjkIzMwyzkFgZpZxiYJA0iRJKyW1SJpeoL1WUrOkJkmNks5IOtbMzNLVbRBIKgHuBM4DyoFLJJV36vYH4MSIqAS+DPyiB2PNzCxFSe5ZXA205O4/jKT5QC2wYluHiNiU138wEEnHWv+pr6/n2tt+knYZ+5TZN19DXV1d2mWY9akkS0MjgDfy9ltzx3Yg6WJJrwAP03FWkHhsbnxdblmpcd26dUlqt27MmzePLWtfS7uMfcaWta8xb968tMsw63NJzghU4FjsdCBiIbBQ0gTgVuCcpGNz4+uBeoCqqqqCfaznBg0fzRGX3p52GfuEt+b5LS7bNyU5I2gFRubtlwFruuocEUuAYyQd2tOxZma25yUJggZgjKTRkgYBU4BF+R0k/Ysk5bbHA4OA9UnGmplZurpdGoqIdklXA48AJcDciFguaVqufQ7wOeBySR8CbcDkiAig4Nh+ei1mZtYLSd4jICIWA4s7HZuTtz0TmJl0rJmZFQ9/s9jMLOMcBGZmGecgMDPLOAeBmVnGOQjMzDLOQWBmlnEOAjOzjHMQmJllnIPAzCzjHARmZhnnIDAzyzgHgZlZxjkIzMwyzkFgZpZxDgIzs4xzEJiZZVyiIJA0SdJKSS2SdrqDt6TLJDXnfp6TdGJe22pJL0pqktTYl8Wbmdnu6/YOZZJKgDuBc+m4GX2DpEURsSKv22vAmRHxD0nnAfXAyXntZ0XE231Yt5mZ9ZEkZwTVQEtErIqILcB8oDa/Q0Q8FxH/yO2+AJT1bZlmZtZfkgTBCOCNvP3W3LGu/Dvwu7z9AB6VtExSXc9LNDOz/pTk5vUqcCwKdpTOoiMIzsg7fHpErJE0HHhM0isRsaTA2DqgDuCoo45KUJaZmfWFJGcErcDIvP0yYE3nTpIqgF8AtRGxftvxiFiT+3MtsJCOpaadRER9RFRFRNVhhx2W/BWYmdluSRIEDcAYSaMlDQKmAIvyO0g6CngA+GJEvJp3fLCkIdu2gYnAS31VvJmZ7b5ul4Yiol3S1cAjQAkwNyKWS5qWa58DzACGAT+TBNAeEVXA4cDC3LH9gHkR8ft+eSVmZtYrSd4jICIWA4s7HZuTt30VcFWBcauAEzsfNzOz4uFvFpuZZZyDwMws4xwEZmYZ5yAwM8s4B4GZWcY5CMzMMs5BYGaWcQ4CM7OMcxCYmWWcg8DMLOMcBGZmGecgMDPLOAeBmVnGOQjMzDLOQWBmlnEOAjOzjHMQmJllXKIgkDRJ0kpJLZKmF2i/TFJz7uc5SScmHWtmZunqNggklQB3AucB5cAlkso7dXsNODMiKoBbgfoejDUzsxQluWdxNdCSu/8wkuYDtcCKbR0i4rm8/i8AZUnHmu1NmpqaqKmpSbuMfUJTUxMMHZl2GUaypaERwBt5+625Y135d+B3PR0rqU5So6TGdevWJSjLbM8aXF5DZWVl2mXsMyorKxlcXpN2GUayMwIVOBYFO0pn0REEZ/R0bETUk1tSqqqqKtjHLE1DKiexmklpl7FPGZJ2AQYkC4JWIP/8rQxY07mTpArgF8B5EbG+J2PNzCw9SZaGGoAxkkZLGgRMARbld5B0FPAA8MWIeLUnY83MLF3dnhFERLukq4FHgBJgbkQslzQt1z4HmAEMA34mCaA9Iqq6GttPr8XMzHohydIQEbEYWNzp2Jy87auAq5KONTOz4uFvFpuZZZyDwMws4xwEZmYZ5yAwM8s4B4GZWcY5CMzMMs5BYGaWcQ4CM7OMcxCYmWWcg8DMLOMcBGZmGecgMDPLOAeBmVnGOQjMzDLOQWBmlnEOAjOzjHMQmJllXKIgkDRJ0kpJLZKmF2g/TtLzkj6Q9I1ObaslvSipSVJjXxVuZmZ9o9tbVUoqAe4EzgVagQZJiyJiRV63d4BrgX/t4mHOioi3d7NWMzPrB0nOCKqBlohYFRFbgPlAbX6HiFgbEQ3Ah/1Qo5mZ9aMkQTACeCNvvzV3LKkAHpW0TFJdV50k1UlqlNS4bt26Hjy8mZntjiRBoALHogfPcXpEjAfOA74maUKhThFRHxFVEVF12GGH9eDhzcxsdyQJglZgZN5+GbAm6RNExJrcn2uBhXQsNZmZWZFIEgQNwBhJoyUNAqYAi5I8uKTBkoZs2wYmAi/1tlgzM+t73X5qKCLaJV0NPAKUAHMjYrmkabn2OZKOABqBocBHkq4HyoFDgYWStj3XvIj4fb+8EjMz65VugwAgIhYDizsdm5O3/RYdS0advQecuDsFmplZ//I3i83MMs5BYGaWcQ4CM7OMcxCYmWWcg8DMLOMcBGZmGecgMDPLOAeBmVnGOQjMzDLOQWBmlnEOAjOzjHMQmJllnIPAzCzjHARmZhnnIDAzyzgHgZlZxjkIzMwyLlEQSJokaaWkFknTC7QfJ+l5SR9I+kZPxpqZWbq6DQJJJcCdwHl03If4Eknlnbq9A1wLzOrFWDMzS1GSexZXAy0RsQpA0nygFlixrUNErAXWSvpsT8cWsnLlSmpqapK+ButCU1MTDB2ZdhlmVuSSBMEI4I28/Vbg5ISPn3ispDqgDoCSgbywan3Cp7AuDR3J4PKatKswsyKXJAhU4FgkfPzEYyOiHqgH2P/IMXHEpbcnfAozM9sdSd4sbgXy1xfKgDUJH393xpqZ2R6QJAgagDGSRksaBEwBFiV8/N0Za2Zme0C3S0MR0S7pauARoASYGxHLJU3Ltc+RdATQCAwFPpJ0PVAeEe8VGttPr8XMzHohyXsERMRiYHGnY3Pytt+iY9kn0VgzMyse/maxmVnGOQjMzDLOQWBmlnEOAjOzjHMQmJllnIPAzCzjHARmZhnnIDAzyzgHgZlZxjkIzMwyzkFgZpZxDgIzs4xzEJiZZZyDwMws4xwEZmYZ5yAwM8s4B4GZWcYlCgJJkyStlNQiaXqBdkmanWtvljQ+r221pBclNUlq7Mvizcxs93V7q0pJJcCdwLlAK9AgaVFErMjrdh4wJvdzMvAfuT+3OSsi3u6zqs3MrM8kOSOoBloiYlVEbAHmA7Wd+tQC/xkdXgAOlnRkH9dqZmb9IEkQjADeyNtvzR1L2ieARyUtk1TX1ZNIqpPUKKlx6/sbEpRlZmZ9odulIUAFjkUP+pweEWskDQcek/RKRCzZqXNEPVAPsP+RYzo/vpmZ9ZMkZwStwMi8/TJgTdI+EbHtz7XAQjqWmszMrEgkCYIGYIyk0ZIGAVOARZ36LAIuz3166BRgQ0S8KWmwpCEAkgYDE4GX+rB+MzPbTd0uDUVEu6SrgUeAEmBuRCyXNC3XPgdYDJwPtADvA1fmhh8OLJS07bnmRcTv+/xVmJlZryV5j4CIWEzHL/v8Y3PytgP4WoFxq4ATd7NGMzPrR/5msZlZxjkIzMwyzkFgZpZxDgIzs4xzEJiZZZyDwMws4xwEZmYZ5yAwM8s4B4GZWcY5CMzMMs5BYGaWcQ4CM7OMcxCYmWWcg8DMLOMcBGZmGecgMDPLOAeBmVnGJQoCSZMkrZTUIml6gXZJmp1rb5Y0PulYMzNLV7dBIKkEuBM4DygHLpFU3qnbecCY3E8d8B89GGtmZilKckZQDbRExKqI2ALMB2o79akF/jM6vAAcLOnIhGPNzCxFSW5ePwJ4I2+/FTg5QZ8RCccCIKmOjrMJgA9en3nBSwlqS9OhwNtpF5GA6+xbrrNvuc6+88neDkwSBCpwLBL2STK242BEPVAPIKkxIqoS1JaavaFGcJ19zXX2LdfZdyQ19nZskiBoBUbm7ZcBaxL2GZRgrJmZpSjJewQNwBhJoyUNAqYAizr1WQRcnvv00CnAhoh4M+FYMzNLUbdnBBHRLulq4BGgBJgbEcslTcu1zwEWA+cDLcD7wJW7GpugrvrevJg9bG+oEVxnX3Odfct19p1e16iIgkv2ZmaWEf5msZlZxjkIzMwyLrUgSHDZihpJGyQ15X5mpFTnXElrJRX8XsOuLq9RRDUWy1yOlPSkpJclLZd0XYE+xTCfSepMfU4llUpaKul/cnV+t0CfVOczYY2pz2VeLSWS/iTpoQJtqf/bzKtlV3X2fD4jYo//0PHG8V+Ao+n4iOn/AOWd+tQAD6VRX6c6JgDjgZe6aD8f+B0d35k4BfjvIqyxWObySGB8bnsI8GqBv/dimM8kdaY+p7k5OjC3PRD4b+CUYprPhDWmPpd5tfxvYF6hetKeyx7U2eP5TOuMYK+59ERELAHe2UWXri6vscckqLEoRMSbEfHH3PZG4GU6vn2erxjmM0mdqcvN0abc7sDcT+dPf6Q6nwlrLAqSyoDPAr/ookvq/zYhUZ09llYQdHVJis5OzZ1S/k7SCXumtB5L+lrSVlRzKWkU8L/o+B9ivqKaz13UCUUwp7klgiZgLfBYRBTdfCaoEYpgLoE7gP8DfNRFe+pzmXMHu64TejifaQVBkktP/BH4REScCPwEeLC/i+qlxJfRSFFRzaWkA4H7gesj4r3OzQWGpDKf3dRZFHMaEVsjopKOb+1XSxrbqUvq85mgxtTnUtIFwNqIWLarbgWO7dG5TFhnj+czrSDo9rIVEfHetlPKiFgMDJR06J4rMbEkl+BIVTHNpaSBdPxyvTciHijQpSjms7s6i2lOczW8CzwFTOrUVBTzCV3XWCRzeTpwkaTVdCxVf1rSrzv1KYa57LbO3sxnWkHQ7aUnJB0hSbntajpqXb/HK+1eV5fXKBrFMpe5Gn4JvBwRP+yiW+rzmaTOYphTSYdJOji3/THgHOCVTt1Snc8kNRbDXEbEjRFRFhGj6Ph99EREfKFTt9T/bSapszfzmeSic30ukl224t+Ar0hqB9qAKZF7S3xPkvQbOt6FP1RSK/BtOt7w2lZnwctrFFmNRTGXdPxv5ovAi7k1Y4CbgKPyak19PklWZzHM6ZHAPeq4AdQA4L8i4iEluPxLkdVYDHNZUJHNZZd2dz59iQkzs4zzN4vNzDLOQWBmlnEOAjOzjHMQmJllnIPAzCzjHARmZhnnIDAzy7j/BwfgnLyp7iVKAAAAAElFTkSuQmCC\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAWsklEQVR4nO3de5hWZb3/8feX4aCgiImVCQIZSngoCUU8bXTnFtgq2sbC+Nlua9s0MS2PmT8zzONWcysoYrLN/cvTdkthomhoaaIEEiqoJJLIRGUagpwZuH9/zIjDODjPyJp5Hr3fr+uai3W413q+z7puPizW4Z5IKSFJykubchcgSWp9hr8kZcjwl6QMGf6SlCHDX5Iy1LZcH9y1a9fUs2fPcn28JH0oPfPMM2+klHbc0v2ULfx79uzJzJkzy/XxkvShFBELi9iPl30kKUOGvyRlyPCXpAwZ/pKUIcNfkjJk+EtShkoK/4gYHBHzImJ+RJzfyPpBEbE0ImbX/VxUfKmSpKI0+Zx/RFQBY4HDgWpgRkRMSim90KDpEymlI1ugRklSwUp5yWs/YH5KaQFARNwFDAMahn+zzJs3j0GDBm3JLiRJH1Ap4b8zsKjefDUwoJF2AyPiWWAxcHZKaW7DBhFxMnAyAFXteHrBm80uWJK05UoJ/2hkWcNf/zUL6JFSWh4RQ4GfA73fs1FK44HxAB126p0++dUrmletJGVu4ZXFXF0v5YZvNdC93nw3as/uN0opLUspLa+bngy0i4iuhVQoSSpcKeE/A+gdEb0ioj0wAphUv0FEfDIiom56v7r9ek1HkipUk5d9Uko1ETEKmAJUARNSSnMj4pS69eOA4cCpEVEDrAJGJH8zvCRVrJKGdK67lDO5wbJx9abHAGOKLU2S1FJ8w1eSMmT4S1KGDH9JypDhL0kZMvwlKUOGvyRlyPCXpAwZ/pKUIcNfkjJk+EtShgx/ScqQ4S9JGTL8JSlDhr8kZcjwl6QMGf6SlCHDX5IyZPhLUoYMf0nKkOEvSRky/CUpQ4a/JGXI8JekDBn+kpQhw1+SMmT4S1KGDH9JypDhL0kZMvwlKUOGvyRlyPCXpAwZ/pKUoZLCPyIGR8S8iJgfEee/T7t9I2J9RAwvrkRJUtGaDP+IqALGAkOAvsDxEdF3M+2uBKYUXaQkqVilnPnvB8xPKS1IKa0F7gKGNdLudOB/gdcLrE+S1AJKCf+dgUX15qvrlm0UETsDxwLj3m9HEXFyRMyMiJnrVy5tbq2SpIKUEv7RyLLUYP464LyU0vr321FKaXxKqX9KqX9Vx+1KLFGSVLS2JbSpBrrXm+8GLG7Qpj9wV0QAdAWGRkRNSunnRRQpSSpWKeE/A+gdEb2APwEjgK/Wb5BS6vXOdETcBvzS4JekytVk+KeUaiJiFLVP8VQBE1JKcyPilLr173udX5JUeUo58yelNBmY3GBZo6GfUvr6lpclSWpJvuErSRky/CUpQ4a/JGXI8JekDBn+kpQhw1+SMmT4S1KGDH9JypDhL0kZMvwlKUOGvyRlyPCXpAwZ/pKUIcNfkjJk+EtShgx/ScqQ4S9JGTL8JSlDhr8kZcjwl6QMGf6SlCHDX5IyZPhLUoYMf0nKkOEvSRky/CUpQ4a/JGXI8JekDBn+kpQhw1+SMmT4S1KGDH9JypDhL0kZKin8I2JwRMyLiPkRcX4j64dFxHMRMTsiZkbEQcWXKkkqStumGkREFTAWOByoBmZExKSU0gv1mk0FJqWUUkTsDdwD9GmJgiVJW66UM//9gPkppQUppbXAXcCw+g1SSstTSqluthOQkCRVrFLCf2dgUb356rplm4iIYyPiJeAB4MTGdhQRJ9ddFpq5fuXSD1KvJKkApYR/NLLsPWf2KaWJKaU+wDHAJY3tKKU0PqXUP6XUv6rjds0qVJJUnFLCvxroXm++G7B4c41TSo8Du0ZE1y2sTZLUQkoJ/xlA74joFRHtgRHApPoNIuIzERF10/2A9sCbRRcrSSpGk0/7pJRqImIUMAWoAiaklOZGxCl168cB/wJ8LSLWAauAr9S7ASxJqjBNhj9ASmkyMLnBsnH1pq8Eriy2NElSS/ENX0nKkOEvSRky/CUpQ4a/JGXI8JekDBn+kpQhw1+SMmT4S1KGDH9JypDhL0kZMvwlKUOGvyRlyPCXpAwZ/pKUIcNfkjJk+EtShgx/ScqQ4S9JGTL8JSlDhr8kZcjwl6QMGf6SlCHDX5IyZPhLUoYMf0nKkOEvSRky/CUpQ4a/JGXI8JekDBn+kpQhw1+SMmT4S1KGSgr/iBgcEfMiYn5EnN/I+pER8Vzdz7SI+FzxpUqSitJk+EdEFTAWGAL0BY6PiL4Nmv0R+IeU0t7AJcD4oguVJBWnlDP//YD5KaUFKaW1wF3AsPoNUkrTUkpL6mafBroVW6YkqUilhP/OwKJ689V1yzbnJODBxlZExMkRMTMiZq5fubT0KiVJhWpbQptoZFlqtGHEodSG/0GNrU8pjafuklCHnXo3ug9JUssrJfyrge715rsBixs2ioi9gZ8AQ1JKbxZTniSpJZRy2WcG0DsiekVEe2AEMKl+g4jYBbgPOCGl9Ifiy5QkFanJM/+UUk1EjAKmAFXAhJTS3Ig4pW79OOAiYAfgxogAqEkp9W+5siVJW6KUyz6klCYDkxssG1dv+hvAN4otTZLUUnzDV5IyZPhLUoYMf0nKkOEvSRky/CUpQ4a/JGXI8JekDBn+kpQhw1+SMmT4S1KGDH9JypDhL0kZMvwlKUOGvyRlyPCXpAwZ/pKUIcNfkjJk+EtShgx/ScqQ4S9JGTL8JSlDhr8kZcjwl6QMGf6SlCHDX5IyZPhLUoYMf0nKkOEvSRky/CUpQ4a/JGXI8JekDBn+kpShksI/IgZHxLyImB8R5zeyvk9EPBURayLi7OLLlCQVqW1TDSKiChgLHA5UAzMiYlJK6YV6zf4OfBs4piWKlCQVq5Qz//2A+SmlBSmltcBdwLD6DVJKr6eUZgDrWqBGSVLBSgn/nYFF9ear65Y1W0ScHBEzI2Lm+pVLP8guJEkFKCX8o5Fl6YN8WEppfEqpf0qpf1XH7T7ILiRJBSgl/KuB7vXmuwGLW6YcSVJrKCX8ZwC9I6JXRLQHRgCTWrYsSVJLavJpn5RSTUSMAqYAVcCElNLciDilbv24iPgkMBPoDGyIiDOBvimlZc0ppnOHNpw+YHt6dGlHNHq1Sc2VSCx8ax03TF/CsjUbyl2OpArRZPgDpJQmA5MbLBtXb/ov1F4O2iKnD9iefrt+irYdtyXC8C9CSokddljG6cClj79Z7nIkVYiKesO3R5d2Bn/BIoK2HTvTo0u7cpciqYJUVPgHYfC3gIjwMpqkTVRU+EuSWofhXwFemvs8Tzz6cLO3O+m4I5n77O9boCJJH3WGfwWYN/d5nnj0kXKXISkjJT3tUw5XXfw95s19vtB97r7HXpx78eVNtvvlfXdzx4Tx1Kxby577fIFjvvJ/GH3uGfzs/qms37CekUd9kavG3spbS/7OjVdfxnbbf4yFC+bTb8BAvn/pNbRp04Zpv3mUm669grVr19C9Ry9GXzOGjp22Yc7sWVx18fmsWrmSdu07cPMd93HjNZexZvVqZs94mhNP+w6HfPEIrvi/5zH/pReoWV/Dqd85n0OPGMrqVau46KzTWPDyPHp9ZndWr15d6PGRlI+KDf9yWfDyPKbcP5GfTnyIdu3acekFZ7HwlfkMOnwIY/7jUtasXsWRxx5H7z59mfHUb5nz7CwmTn2anbp151snDGfqg/fTf+BB3HL91dx850Q6duzEhBuv4/ZbbuSkb53JuaedyFVjJ7Dn5/ux/O1lbLV1R7511gXMfe73XPCj/wDg+itGs9+BBzP6mjEsW7qUkUf9IwMO/gfu/X+3sdXWHbn3kSf5w4tzGDFkUHkPlqQPrYoN/1LO0FvC9Cd/w4vPPcvIIw8DYPXq1Xys645888xz+eqRh9G+w1acN/rKje33/Fw/uvXoCcDgo/+F3894mvYdOrDg5Xl8/djBAKxbt469++3Lq6+8zI4f/wR7fr4fANts27nRGp56/DF+/ciD3H7zGADWrlnNX/5Uzazp0zj+xG8CsNtn96T3Z/dokWMg6aOvYsO/XFKCo44bwRnn/2CT5W+8/ldWrlhBTU0Na9aspmPHTgDveTQ1CFJK7H/wIK4ce+sm6/7w4hwo4VHWROLa8bfTc9fe71nno7CSiuAN3wYGHHgIv3pgEm++8TcAli5ZwuLq1xh93pmcdvYFDD1mONdddvHG9nNmz6L6tYVs2LCBKfdPZJ/99mfvfvsye+Z0XvvjAgBWrVrJqwvm02vX3fjbX//CnNmzAFix/G1qamrotM02rFyxfOM+DzjkMO74r/GkVDt46otzngOg34ADmDzxfwB4+aUXePnFuS1+PCR9NHnm38Cuu/XhtHO+z6kjv8SGDRto264dg/5pKFVt2zL02ONYv349XzvmCKY/+Tht2rRh7y/sy39e/kPmz3uBfgMGctjgI2nTpg2jr72R80d9g7Vr1wAw6pzv0/PTn+GqsRO44qLzWLN6FR222prxd05k34EHM2HsdXz5iIM58bTvcPIZ53DVD7/H8MMPJKXEp7rvwpjb7ubLJ5zIRWedxvDDD2T3PfbaePlIkpor3jm7bG0dduqddvrX6zZZdsvRO/GJXT5dlno+iBlP/Zaf3nwDY267u9ylNOmvry3g3yf9udxlSNpCC6888pmUUv8t3Y+XfSQpQ1722QL7DjyIfQceVO4yJKnZPPOXpAwZ/pKUIcNfkjJk+EtShir6hu/RY54sdH+TRh3YZJs/LXqN07/+Fe6b+tQmy8defRlfGHAA+x88qNHtHn3oAXp8eld23a1PEaVKUovyzL9Ep519wWaDH+CxKQ+w4OV5hXxWTU1NIfuRpM2p6DP/ctmwYQM/PPcMZj/zOz7+iZ34z1t/xqUXnMUhXzyCw/95GNddfjG/eeQhqqqqGHjIYfzjkCP59SMPMnP6k9xy/dVcc/PtrFixnB9977usXrWSbj16MfrqMXTu0oU5s2dx8Tmns3XHTuyz7wB++9ivuG/qU/zinjt44tGHWbNmNatWruT6CXdwxkkjWbb0LWrWrWPUORdy6BFD+dOi1/jWCcPZZ9/9eW7WDHbvuyfDvjySm669nL+/8QaXXT+evfb5QrkPoaQK55l/I1774yt85V+/wcSpT9F5u+341YOTNq5bumQJjz70APdNfYp7H3mSf//22Xy+/wAGHT6E735/NPdMeYLuPXtx4ZmncOb3LubeR56kd5++jLuudiTQi846jQsvv5b//sXDtKmq2uRzn31mBj+69iZ+cvck2nfYih/f8t/c/eBv+Mk993PNJRduHOtn0asLGHniN7n3kSf54ysvM/nn93LbfQ/x3Qsv4dYx17begZL0oWX4N2Ln7j3os8deAHx2r8+xeNGijes6bbstHTp04OJzvs2vHryfrbfe+j3bv71sKW8vW0r/gbX3GI4efjzPTJ/GsqVLWbliOZ/vPwCAoccM32S7/Q8exHbbbw9ASonrr7yE4YcfyDePP4bX//Jn3vzb6xvr6/3ZPWjTpg277taHAQcdQkTQu09fFle/VvwBkfSRY/g3ol379hunq9pUsX79u9fg27Zty8/un8oXhx7FY1Me4NQThje2i8Y1MY7S1h07bpyePPF/WPLmm9w5+dfcM+UJdthxR9asWfOe+tpEG9q37wBAtGnj/QJJJTH8m2nliuW8/fYyDj7snzj3B5dv/FWTHbfZhhXLa4dl3rbzdnTerguzpk8Dan8tZP8BB9K5Sxc6dtqG52bNAOChX9y32c9Z/vYyPta1K+3ateN3055gcfWizbaVpOaq6Bu+pTya2dpWLF/OGSeNZO2a1aSUOOcHlwEw+OgvMfq8M7njv27mmnE/5ZIf3/TuDd9dejL6mrEAXHz1DYw+9wy27tiJ/gMPZNvOjf82r6HHHse3/+14jh96KLvvsRe9PrNbq31HSR99DuncylauWE7HTtsAcOvYH/PG63/lvB9e0eKf65DO0kdDUUM6V/SZ/0fR41MfZsLYH1NTU8OnunVn9LU3lrskSRky/FvZ4KO/xOCjv1TuMiRlrqJu+CYS5boM9VGWUiLhcZX0rooK/4VvraNm5TL/AShQSomalctY+Na6cpciqYJU1GWfG6Yv4XSgR5c3CKLc5XwkJBIL31rHDdOXlLsUSRWkosJ/2ZoNXPr4m+UuQ5I+8kq67BMRgyNiXkTMj4jzG1kfEXF93frnIqJf8aVKkorSZPhHRBUwFhgC9AWOj4i+DZoNAXrX/ZwM3FRwnZKkApVy5r8fMD+ltCCltBa4CxjWoM0w4PZU62mgS0TsVHCtkqSClHLNf2eg/sAy1cCAEtrsDGzySmlEnEzt/wwA1iy88sg5zaq2PLoCb5S7iBJYZ7E+DHV+GGoE6yza7kXspJTwb+yxm4bPYpbShpTSeGA8QETMLOIV5ZZmncWyzuJ8GGoE6yxaRMwsYj+lXPapBrrXm+8GLP4AbSRJFaKU8J8B9I6IXhHRHhgBTGrQZhLwtbqnfvYHlqaUHEVMkipUk5d9Uko1ETEKmAJUARNSSnMj4pS69eOAycBQYD6wEvi3Ej57/AeuunVZZ7GsszgfhhrBOotWSJ1lG9JZklQ+FTW2jySpdRj+kpShFgn/EoaDGBQRSyNidt3PRaVu24o1nlOvvjkRsT4iPla37tWIeL5uXSGPXb1PnRMi4vWIaPSdiPcbWqO1jmWJdY6sq++5iJgWEZ+rt66Sjmcl9M2maqyUvtk9Ih6LiBcjYm5EnNFIm7L3zxLrLHv/LLHO4vpnSqnQH2pvCr8CfBpoDzwL9G3QZhDwyw+ybWvV2KD9UcCj9eZfBboWXddmPvsQoB8wZzPrhwIPUvuuxf7A9NY8ls2o8wBg+7rpIe/UWYHHs6x9s5QaK6hv7gT0q5veFvhDI3/Xy94/S6yz7P2zxDoL658tceZfynAQLbFtS9Z4PHBnC9TRpJTS48Df36fJ5obWaK1jWVKdKaVpKaV3xpV+mtp3QVpdCcdzc1rteDazxnL2zT+nlGbVTb8NvEjtm/31lb1/llJnJfTPEo/n5jT7eLZE+G9uqIeGBkbEsxHxYETs0cxtW6tGIqIjMBj433qLE/BwRDwTtUNWlNPmvktrHcsP4iRqzwbfUUnHE8rbN0tWSX0zInoC+wDTG6yqqP75PnXWV/b+2USdhfTPlhjPv5ShHmYBPVJKyyNiKPBzakcELWmYiAI053OOAp5MKdU/EzswpbQ4Ij4OPBIRL9WdrZXD5r5Lax3LZomIQ6n9y3VQvcWVdDzL3TeboyL6ZkRsQ+0/QGemlJY1XN3IJmXpn03U+U6bsvfPJuosrH+2xJl/k0M9pJSWpZSW101PBtpFRNdStm2tGusZQYP/VqeUFtf9+Towkdr/cpXL5r5LxQ25ERF7Az8BhqWUNv7Wnko6nhXQN5uj7H0zItpRG1Q/Synd10iTiuifJdRZEf2zqToL7Z8tcNOiLbAA6MW7Nx72aNDmk7z7gtl+wGvU/svV5LatVWNdu+2ovfbaqd6yTsC29aanAYOLrrFBHT3Z/A3Kf2bTG2q/a853bMU6d6H2DfADGiyvtONZ1r5ZSo2V0jfrjsvtwHXv06bs/bPEOsveP0uss7D+Wfhln1TacBDDgVMjogZYBYxItd+m0W3LVCPAscDDKaUV9Tb/BDAxIqD2gN+RUnqo6BrfERF3UnuHv2tEVAM/ANrVq7PRoTU29x3LWOdFwA7AjXXHribVjqBYacezrH2zxBqhAvomcCBwAvB8RMyuW3YBtUFaSf2zlDoroX+WUmdh/dPhHSQpQ77hK0kZMvwlKUOGvyRlyPCXpAwZ/pKUIcNfkjJk+EtShv4/2z/BfwfSC9YAAAAASUVORK5CYII=\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAX7klEQVR4nO3df3BV5Z3H8fc3IRBBQAulMoSV6GIxQpqlEAFdTMuPRhbNMm2HX10tXWGwxepM6zQ6lrbKVhiZ1uKPpUyLuK0UdhBpBulqqbIUlSXBppEfxqYoGtGiUBAkiMHv/pErc4035iS5l5ub5/Oaycz58dxzv899mA8nzz3nxNwdERHp+rLSXYCIiJwdCnwRkUAo8EVEAqHAFxEJhAJfRCQQ3dL1xv379/chQ4ak6+1FRDLSzp0733b3T7fntWkL/CFDhlBVVZWutxcRyUhmtr+9r9WUjohIIBT4IiKBUOCLiAQibXP4ItI1vP/++9TX13Py5Ml0l9Kl5ObmkpeXR05OTtKOqcAXkQ6pr6+nd+/eDBkyBDNLdzldgrtz6NAh6uvryc/PT9pxNaUjIh1y8uRJ+vXrp7BPIjOjX79+Sf+tKVLgm1mpmdWaWZ2ZlSfYf6uZVcd+dpnZaTP7VFIrFZFOS2GffKn4TFsNfDPLBh4ArgYKgJlmVhDfxt3vcfcidy8CbgP+190PJ71aERFptyhz+MVAnbvvAzCzNUAZsKeF9jOB37R20NraWkpKSiKWKZ3NrFmzmDdvXrrLEEmb6upqDhw4wJQpU9r0upKSEpYuXcqoUaNSVFnLogT+IOC1uPV64PJEDc2sJ1AKLGhh/zygKSWyc9i+71BbapVO4tTBlwEU+BK06upqqqqq2hz46RQl8BNNJLX0Z7KuAZ5paTrH3VcAKwB6DBzqF8xaHKlI6VzeXP2xr3FE0u7Xv/41y5Yt49SpU1x++eV84xvfYO7cuezYsYPTp09TXFzM2rVrefvtt1m4cCH9+vWjtraW8ePH8+CDD5KVlcWTTz7JD37wA9577z0uvvhiHnroIc4991wqKyu5+eabeffdd+nRowe///3vWbhwIQ0NDWzbto3bbruNqVOnctNNN/HCCy/Q2NjID3/4Q8rKymhoaGDOnDns2bOHSy+9lIaGhrR9RlECvx4YHLeeBxxooe0MIkzniEjXdMstt1BdXZ3UYxYVFXHvvfd+Ypu9e/eydu1annnmGXJycvjmN79JbW0t1157LXfccQcNDQ187WtfY/jw4WzZsoUdO3awZ88eLrzwQkpLS1m/fj0lJSUsWrSIzZs306tXL5YsWcJPfvITysvLmT59OmvXrmX06NG888479OzZkzvvvJOqqiruv/9+AG6//Xa++MUvsnLlSo4cOUJxcTETJ07k5z//OT179qSmpoaamhpGjhyZ1M+nLaIEfiUw1MzygddpCvVZzRuZWV/gKuBrSa1QRKQVf/jDH9i5cyejR48GoKGhgQEDBrBw4UJGjx5Nbm4uy5YtO9O+uLiYiy66CICZM2eybds2cnNz2bNnD1dccQUAp06dYuzYsdTW1jJw4MAzx+7Tp0/CGp588kkqKipYunQp0HS56quvvsrWrVv59re/DUBhYSGFhYWp+RAiaDXw3b3RzBYATwDZwEp3321m82P7l8eaTgOedPd3U1atiHRqrZ2Jp4q7c/3113P33Xd/ZPubb77J8ePHef/99zl58iS9evUCPn7Jo5nh7kyaNInf/OajkxQ1NTWRLpF0dx599FE++9nPfmxfZ7lsNdJ1+O6+yd0vcfeL3f0/YtuWx4U97r7K3WekqlARkZZMmDCBdevWcfDgQQAOHz7M/v37mTdvHnfddRezZ8/me9/73pn2O3bs4OWXX+aDDz5g7dq1XHnllYwZM4ZnnnmGuro6AE6cOMFLL73EsGHDOHDgAJWVlQAcO3aMxsZGevfuzbFjx84c80tf+hL33Xcf7k1fcf7pT38CYPz48TzyyCMA7Nq1i5qamtR/IC3QoxVEJOMVFBSwaNEiJk+ezAcffEBOTg5lZWV069aNWbNmcfr0acaNG8dTTz1FVlYWY8eOpby8nBdeeIHx48czbdo0srKyWLVqFTNnzuS9994DYNGiRVxyySWsXbuWm266iYaGBs455xw2b97MF77wBRYvXkxRURG33XYb3//+97nlllsoLCzE3RkyZAgbN27kxhtvZM6cORQWFlJUVERxcXHaPif78H+js63HwKE+8Pp70/Le0jFvri5nzEX92LJlS7pLkU5g7969XHrppekuI7ItW7awdOlSNm7cmO5SWpXoszWzne7erov49SwdEZFAaEpHRIJSUlIS7F3+OsMXEQmEAl9EJBAKfBGRQCjwRUQCoS9tRSSphpQ/ntTjvbL4X1pv88orTJ06lV27dn1k+8KFCxk/fjwTJ05M+LoNGzZwySWXUFBQkHB/V6MzfBHpsu68884Wwx6aAn/Pnpb+tEfbNDY2JuU4qaTAF5Eu4fTp08ydO5fLLruMyZMn09DQwNe//nXWrVsHQHl5OQUFBRQWFvLd736XZ599loqKCm699VaKior461//SnV1NWPGjKGwsJBp06bx97//HYDKykoKCwsZO3Yst956K8OHDwdg1apVfPWrX+Waa65h8uTJHD9+nAkTJjBy5EhGjBjBb3/7W6DpN5Bhw4Zxww03MHz4cGbPns3mzZu54oorGDp0KDt27Dgrn5ECX0S6hL/85S9861vfYvfu3Zx33nk8+uijZ/YdPnyYxx57jN27d1NTU8Mdd9zBuHHjuPbaa7nnnnuorq7m4osv5rrrrmPJkiXU1NQwYsQIfvSjHwEwZ84cli9fznPPPUd2dvZH3ve5557j4Ycf5qmnniI3N5fHHnuM559/nqeffprvfOc7Z56tU1dXx80330xNTQ0vvvgiq1evZtu2bSxdupQf//jHZ+UzUuCLSJeQn59PUVERAJ///Od55ZVXzuzr06cPubm53HDDDaxfv56ePXt+7PVHjx7lyJEjXHXVVQBcf/31bN26lSNHjnDs2DHGjRsHNP15z3iTJk3iU5/6FND0xMzbb7+dwsJCJk6cyOuvv87f/va3M/WNGDGCrKwsLrvsMiZMmICZMWLEiI/UmkoKfBHpEnr06HFmOTs7+yNz6t26dWPHjh18+ctfZsOGDZSWlkY+bmvPG/vwkcsAjzzyCG+99RY7d+6kurqaz3zmM5w8efJj9WVlZZ1Zz8rKOmvz/wp8Eenyjh8/ztGjR5kyZQr33nvvmb/KFf+I4759+3L++efzxz/+EYBf/epXXHXVVZx//vn07t2b7du3A7BmzZoW3+fo0aMMGDCAnJwcnn76afbv35/ajrWRLssUkaSKchnl2Xbs2DHKyso4efIk7s5Pf/pTAGbMmMHcuXNZtmwZ69at4+GHH2b+/PmcOHGCiy66iIceegiAX/7yl8ydO5devXpRUlJC3759E77P7Nmzueaaaxg1ahRFRUUMGzbsrPUxCj0eWdpMj0eWeJn2eOT2OH78OOeeey4Aixcv5o033uBnP/tZyt832Y9H1hm+iEgrHn/8ce6++24aGxu58MILWbVqVbpLahcFvohIK6ZPn8706dPTXUaH6UtbEemwdE0Nd2Wp+EwjBb6ZlZpZrZnVmVl5C21KzKzazHab2f8mt0wR6axyc3M5dOiQQj+J3J1Dhw6Rm5ub1OO2OqVjZtnAA8AkoB6oNLMKd98T1+Y84EGg1N1fNbMBSa1SRDqtvLw86uvreeutt9JdSpeSm5tLXl5eUo8ZZQ6/GKhz930AZrYGKAPinzg0C1jv7q8CuPvBpFYpIp1WTk4O+fn56S5DIogypTMIeC1uvT62Ld4lwPlmtsXMdprZdYkOZGbzzKzKzKpOnzjavopFRKRdopzhW4JtzSfrugGfByYA5wDPmdl2d3/pIy9yXwGsgKbr8NteroiItFeUwK8HBset5wEHErR5293fBd41s63A54CXEBGRTiHKlE4lMNTM8s2sOzADqGjW5rfAP5tZNzPrCVwO7E1uqSIi0hGtnuG7e6OZLQCeALKBle6+28zmx/Yvd/e9ZvY/QA3wAfALd9/V8lFFRORsi3SnrbtvAjY127a82fo9wD3JK01ERJJJd9qKiARCgS8iEggFvohIIBT4IiKBUOCLiARCgS8iEggFvohIIBT4IiKBUOCLiARCgS8iEggFvohIIBT4IiKBUOCLiARCgS8iEggFvohIIBT4IiKBUOCLiARCgS8iEggFvohIIBT4IiKBUOCLiAQiUuCbWamZ1ZpZnZmVJ9hfYmZHzaw69rMw+aWKiEhHdGutgZllAw8Ak4B6oNLMKtx9T7Omf3T3qSmoUUREkiDKGX4xUOfu+9z9FLAGKEttWSIikmxRAn8Q8Frcen1sW3NjzezPZvY7M7ss0YHMbJ6ZVZlZ1ekTR9tRroiItFeUwLcE27zZ+vPAhe7+OeA+YEOiA7n7Cncf5e6jsnv2bVOhIiLSMVECvx4YHLeeBxyIb+Du77j78djyJiDHzPonrUoREemwKIFfCQw1s3wz6w7MACriG5jZBWZmseXi2HEPJbtYERFpv1av0nH3RjNbADwBZAMr3X23mc2P7V8OfAW40cwagQZghrs3n/YREZE0ajXw4cw0zaZm25bHLd8P3J/c0kREJJl0p62ISCAU+CIigVDgi4gEQoEvIhIIBb6ISCAU+CIigVDgi4gEQoEvIhIIBb6ISCAU+CIigVDgi4gEQoEvIhIIBb6ISCAU+CIigVDgi4gEQoEvIhIIBb6ISCAU+CIigVDgi4gEQoEvIhKISIFvZqVmVmtmdWZW/gntRpvZaTP7SvJKFBGRZGg18M0sG3gAuBooAGaaWUEL7ZYATyS7SBER6bgoZ/jFQJ2773P3U8AaoCxBu5uAR4GDSaxPRESSJErgDwJei1uvj207w8wGAdOA5Z90IDObZ2ZVZlZ1+sTRttYqIiIdECXwLcE2b7Z+L/A9dz/9SQdy9xXuPsrdR2X37BuxRBERSYZuEdrUA4Pj1vOAA83ajALWmBlAf2CKmTW6+4ZkFCkiIh0XJfArgaFmlg+8DswAZsU3cPf8D5fNbBWwUWEvItK5tBr47t5oZgtouvomG1jp7rvNbH5s/yfO24uISOcQ5Qwfd98EbGq2LWHQu/vXO16WiIgkm+60FREJhAJfRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAqHAFxEJRKTAN7NSM6s1szozK0+wv8zMasys2syqzOzK5JcqIiId0a21BmaWDTwATALqgUozq3D3PXHN/gBUuLubWSHw38CwVBQsIiLtE+UMvxioc/d97n4KWAOUxTdw9+Pu7rHVXoAjIiKdSpTAHwS8FrdeH9v2EWY2zcxeBB4HvpHoQGY2LzblU3X6xNH21CsiIu0UJfAtwbaPncG7+2PuPgz4V+CuRAdy9xXuPsrdR2X37NumQkVEpGOiBH49MDhuPQ840FJjd98KXGxm/TtYm4iIJFGUwK8EhppZvpl1B2YAFfENzOwfzcxiyyOB7sChZBcrIiLt1+pVOu7eaGYLgCeAbGClu+82s/mx/cuBLwPXmdn7QAMwPe5LXBER6QRaDXwAd98EbGq2bXnc8hJgSXJLExGRZNKdtiIigVDgi4gEQoEvIhIIBb6ISCAU+CIigVDgi4gEQoEvIhIIBb6ISCAU+CIigVDgi4gEQoEvIhIIBb6ISCAU+CIigVDgi4gEQoEvIhIIBb6ISCAU+CIigVDgi4gEQoEvIhIIBb6ISCAiBb6ZlZpZrZnVmVl5gv2zzawm9vOsmX0u+aWKiEhHtBr4ZpYNPABcDRQAM82soFmzl4Gr3L0QuAtYkexCRUSkY6Kc4RcDde6+z91PAWuAsvgG7v6su/89trodyEtumSIi0lHdIrQZBLwWt14PXP4J7f8d+F1HihKR1FixYgWrV69OdxmSJlEC3xJs84QNzb5AU+Bf2cL+ecA8gOw+n45Yoogky+rVq9m6vYruA/LTXYqkQZTArwcGx63nAQeaNzKzQuAXwNXufijRgdx9BbH5/R4Dhyb8T0NEUqv7gHwumLU43WVIO+1fMrXdr40yh18JDDWzfDPrDswAKuIbmNk/AOuBf3P3l9pdjYiIpEyrZ/ju3mhmC4AngGxgpbvvNrP5sf3LgYVAP+BBMwNodPdRqStbRETaKsqUDu6+CdjUbNvyuOUbgBuSW5qIiCST7rQVEQmEAl9EJBAKfBGRQCjwRUQCocAXEQmEAl9EJBAKfBGRQCjwRUQCocAXEQmEAl9EJBAKfBGRQCjwRUQCocAXEQmEAl9EJBAKfBGRQCjwRUQCocAXEQmEAl9EJBAKfBGRQCjwRUQCocAXEQlEpMA3s1IzqzWzOjMrT7B/mJk9Z2bvmdl3k1+miIh0VLfWGphZNvAAMAmoByrNrMLd98Q1Owx8G/jXVBQpIiId12rgA8VAnbvvAzCzNUAZcCbw3f0gcNDM/iUlVUqnU11dTUlJSbrLkDaqrq6GPoPTXYakSZTAHwS8FrdeD1zenjczs3nAPIDsPp9uzyGkE+hVUMK7e7awfd+hdJcibdVnML0KStJdhaRJlMC3BNu8PW/m7iuAFQA9Bg5t1zEk/XoXldK7qDTdZYhIG0X50rYeiP8dMA84kJpyREQkVaIEfiUw1Mzyzaw7MAOoSG1ZIiKSbK1O6bh7o5ktAJ4AsoGV7r7bzObH9i83swuAKqAP8IGZ3QIUuPs7qStdRETaIsocPu6+CdjUbNvyuOU3aZrqERGRTkp32oqIBEKBLyISCAW+iEggFPgiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBCJS4JtZqZnVmlmdmZUn2G9mtiy2v8bMRia/VBER6YhWA9/MsoEHgKuBAmCmmRU0a3Y1MDT2Mw/4zyTXKSIiHRTlDL8YqHP3fe5+ClgDlDVrUwb8lzfZDpxnZgOTXKuIiHRAtwhtBgGvxa3XA5dHaDMIeCO+kZnNo+k3AID39i+ZuqtN1WaW/sDb6S4ihdS/zNWV+wZdv3+fbe8LowS+Jdjm7WiDu68AVgCYWZW7j4rw/hlJ/ctsXbl/XblvEEb/2vvaKFM69cDguPU84EA72oiISBpFCfxKYKiZ5ZtZd2AGUNGsTQVwXexqnTHAUXd/o/mBREQkfVqd0nH3RjNbADwBZAMr3X23mc2P7V8ObAKmAHXACWBOhPde0e6qM4P6l9m6cv+6ct9A/WuRuX9sql1ERLog3WkrIhIIBb6ISCBSHvgRHstQYmZHzaw69rMw1TUli5mtNLODZpbwfoJMf+REhP5l8tgNNrOnzWyvme02s5sTtMnY8YvYv0wev1wz22Fmf47170cJ2mTy+EXpX9vHz91T9kPTl7x/BS4CugN/BgqatSkBNqayjhT2bzwwEtjVwv4pwO9ouk9hDPB/6a45yf3L5LEbCIyMLfcGXkrwbzNjxy9i/zJ5/Aw4N7acA/wfMKYLjV+U/rV5/FJ9hh/lsQwZy923Aoc/oUlGP3IiQv8ylru/4e7Px5aPAXtpujs8XsaOX8T+ZazYmByPrebEfppfgZLJ4xelf22W6sBv6ZELzY2N/eryOzO7LMU1nU1R+5/JMn7szGwI8E80nUXF6xLj9wn9gwwePzPLNrNq4CDwe3fvUuMXoX/QxvFLdeBHeeTC88CF7v454D5gQ4prOpsiPXIig2X82JnZucCjwC3u/k7z3QleklHj10r/Mnr83P20uxfRdGd/sZkNb9Yko8cvQv/aPH6pDvxWH7ng7u98+KuLu28Ccsysf4rrOlu69CMnMn3szCyHpjB8xN3XJ2iS0ePXWv8yffw+5O5HgC1AabNdGT1+H2qpf+0Zv1QHfquPZTCzC8zMYsvFsZoOpbius6VLP3Iik8cuVvcvgb3u/pMWmmXs+EXpX4aP36fN7LzY8jnARODFZs0yefxa7V97xi/K0zLbzaM9luErwI1m1gg0ADM89hV0Z2dmv6Hpm/L+ZlYP/ICmL1c+7Ft7HjnRaUToX8aOHXAF8G/AC7F5UoDbgX+ALjF+UfqXyeM3EHjYmv5AUxbw3+6+0Tr+yJfOIkr/2jx+erSCiEggdKetiEggFPgiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBOL/ATZnvW/lEAPSAAAAAElFTkSuQmCC\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"def f_inv_discrete(prob_list,p):\n",
|
|
" assert np.isclose(np.sum(prob_list), 1), \"The probabilities should sum to one.\"\n",
|
|
" \n",
|
|
" p_cum = 0\n",
|
|
" i = 0\n",
|
|
" while p_cum < p:\n",
|
|
" p_cum += prob_list[i]\n",
|
|
" i += 1\n",
|
|
" return i\n",
|
|
"\n",
|
|
"# plotting\n",
|
|
"f_X = lambda prob_list, x: prob_list[np.rint(x).astype(int) - 1] if np.rint(x) in range(1, len(prob_list) + 1) else 0\n",
|
|
"\n",
|
|
"# Input parameters as list for flexibility in testing.\n",
|
|
"for prob_list in [[.1, .3, .2, .4], [0.5, 0.5], [0.7, 0.1, 0.2]]:\n",
|
|
" alpha, b = params\n",
|
|
" pdf = lambda x: f_X(prob_list, x)\n",
|
|
" samples = [inversion_sample(lambda p: f_inv_discrete(prob_list, p)) for _ in range(100000)]\n",
|
|
" compare_plot(samples, pdf, .5, len(prob_list) + .5, len(prob_list))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 8,
|
|
"id": "3c691f0a",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "b11d87e414ba9dfe2741d73dd95a2f12",
|
|
"grade": true,
|
|
"grade_id": "cell-140af6b31464fbef",
|
|
"locked": true,
|
|
"points": 15,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"assert f_inv_discrete([0.5,0.5],0.4)==1\n",
|
|
"assert f_inv_discrete([0.5,0.5],0.8)==2\n",
|
|
"assert f_inv_discrete([0,0,1],0.1)==3"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "47546d37",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "32dd38f0f963c6132fcbe3ef1f5b9682",
|
|
"grade": false,
|
|
"grade_id": "cell-49fd13dc534dfa28",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"## Central limit theorem? \n",
|
|
"__(35 Points)__\n",
|
|
"\n",
|
|
"In this exercise we will have a closer look at central limits of the Pareto distribution, for which you implemented a random sampler in the previous exercise. By performing the appropriate integrals it is straightforward to show that \n",
|
|
"\n",
|
|
"$$ \n",
|
|
"\\mathbb{E}[X] = \\begin{cases} \\infty & \\text{for }\\alpha \\leq 1 \\\\ \\frac{\\alpha b}{\\alpha - 1} & \\text{for }\\alpha > 1 \\end{cases}, \\qquad \\operatorname{Var}(X) = \\begin{cases} \\infty & \\text{for }\\alpha \\leq 2 \\\\ \\frac{\\alpha b^2}{(\\alpha - 1)^2(\\alpha-2)} & \\text{for }\\alpha > 2 \\end{cases}.\n",
|
|
"$$\n",
|
|
"\n",
|
|
"This shows in particular that the distribution is **heavy tailed**, in the sense that some moments $\\mathbb{E}[X^k]$ diverge."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "ccae582d",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "e6d5659ef88eccfb693b35a088d0d50f",
|
|
"grade": false,
|
|
"grade_id": "cell-a05e255c144ef6c5",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"__(a)__ Write a function `sample_Zn` that produces a random sample for $Z_n= \\frac{\\sqrt{n}}{\\sigma_X}(\\bar{X}_n - \\mathbb{E}[X])$ given $\\alpha>2$, $b>0$ and $n\\geq 1$. Visually verify the central limit theorem for $\\alpha = 4$, $b=1$ and $n=1000$ by comparing a histogram of $Z_n$ to the standard normal distribution (you may use `compare_plot`). __(10 pts)__"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 9,
|
|
"id": "82fe6efd",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "177917ec75361799067d6c23a28569cd",
|
|
"grade": false,
|
|
"grade_id": "cell-b7186322b09717f8",
|
|
"locked": false,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAcbklEQVR4nO3df3CU9bn38fcFJqamQS2SsQNoQqulkaRbGlLAFraiNP7kOE8dcGn99Wikrb9mqo/UY5mpMhpHp1XR89BMa/Vp3cIZ1JZRjnpQo9VKIdhtgCiWJlii0nBixfAjRuT7/LHZZbPkxxJ2c+/e+3nNMGb3vnf34jZcufK9v9/ra845REQk943yOgAREUkPJXQREZ9QQhcR8QkldBERn1BCFxHxiWO8+uCTTjrJlZWVefXxIiI5aePGjf/jnBvX3zHPEnpZWRlNTU1efbyISE4ys3cGOqYhFxERn1BCFxHxCSV0ERGf8GwMXUT84ZNPPqG9vZ3u7m6vQ/GVoqIiJkyYQEFBQcqvUUIXkaPS3t5OSUkJZWVlmJnX4fiCc47Ozk7a29spLy9P+XUachGRo9Ld3c3YsWOVzNPIzBg7duwR/9ajhC4iR03JPP2Gc02V0EVEfEIJXURkGCKRCGvWrDni1wWDwYwtqtRNUZE0KVv8TPzr7fXnexiJjIRIJEJTUxPnnXee16HEqUIXEV/47W9/S01NDYFAgGuvvZY///nPVFVV0d3dzd69eznjjDPYvHkzjY2NzJo1i4svvpiKigoWLVrEwYMHAXj++eeZMWMGU6dO5ZJLLmHPnj0AbNiwgZkzZ/KVr3yFmpoadu/ezZIlS1i5ciWBQICVK1eyd+9errrqKqZNm8ZXv/pV/vCHPwCwf/9+FixYQFVVFfPnz2f//v0Zuwaq0EUkbW666SYikUha3zMQCHD//fcPes6bb77JypUree211ygoKOAHP/gBW7du5aKLLuL2229n//79fPe732XKlCk0Njayfv16WlpaOPXUU6mtreXJJ58kGAyydOlS1q5dS3FxMffccw8/+9nPWLx4MfPnz2flypVMmzaNjz76iOOOO4477riDpqYmHnroIQBuu+02zjrrLB555BE+/PBDampqOPvss/nFL37BcccdR3NzM83NzUydOjWt1yeRErqI5LwXXniBjRs3Mm3aNCBaFZeWlrJkyRKmTZtGUVERDz74YPz8mpoaJk2aBMCll17Kq6++SlFRES0tLZx55pkA9PT0MGPGDLZu3crnP//5+HuPGTOm3xief/55Vq9ezX333QdEp3P+4x//4JVXXuGGG24AoKqqiqqqqsxcBJTQRSSNhqqkM8U5x+WXX87dd9/d5/mdO3eyZ88ePvnkE7q7uykuLgYOnxJoZjjnOOecc/jd737X51hzc3NKUwidczzxxBN86UtfOuzYSE3r1Bi6SD/KFj8T/yPZb86cOaxatYqOjg4APvjgA9555x3q6uq48847WbhwIbfeemv8/PXr19PW1sbBgwdZuXIl3/jGN5g+fTqvvfYa27ZtA2Dfvn28/fbbTJ48mffee48NGzYA0NXVxYEDBygpKaGrqyv+nt/+9rdZtmwZzjkA/vKXvwAwa9YsHn/8cQA2b95Mc3Nzxq6DKnQRyXkVFRUsXbqUuXPncvDgQQoKCpg3bx7HHHMMoVCITz/9lJkzZ/Liiy8yatQoZsyYweLFi9m0aVP8BumoUaN49NFHufTSS/n4448BWLp0KaeffjorV67k+uuvZ//+/XzmM59h7dq1fOtb36K+vp5AIMCPf/xjfvKTn3DTTTdRVVWFc46ysjKefvppvv/973PllVdSVVVFIBCgpqYmY9fBYj9NRlp1dbXTBheSrYYzBTFfpy2++eabfPnLX/Y6jJQ1NjZy33338fTTT3sdypD6u7ZmttE5V93f+RpyERHxCQ25SF5JdxWtMfbcEwwGCQaDXoeREarQRUR8IqWEbma1ZrbVzLaZ2eJBzptmZp+a2XfSF6KIiKRiyIRuZqOBh4FzgQrgUjOrGOC8e4Dn0h2kiIgMLZUKvQbY5pxrdc71ACuAef2cdz3wBNCRxvhERCRFqdwUHQ/sSHjcDnw98QQzGw9cDJwFTBvojcysDqgDOOWUU440VpGMaGhoIBwO93luZ2vnoeOT3qWurm6kw8pZ6b5RnMrN6+3bt3PBBRewefPmPs8vWbKEWbNmcfbZZ/f7ut///vecfvrpVFQcNuiQk1Kp0Ptbs5o8ef1+4Fbn3KeDvZFzrsE5V+2cqx43blyKIYpkVjgcHrChVE9H22HJXnLHHXfcMWAyh2hCb2lpSctnHThwIC3vczRSSejtwMSExxOA95LOqQZWmNl24DvAf5jZv6UjQJFM6Io8G5++FolECAQCNDY2xv+cHKrn5FA9haXlRCKR+LkNDQ1ehy4D+PTTT7nmmms444wzmDt3Lvv37+eKK65g1apVACxevJiKigqqqqq4+eab+dOf/sTq1au55ZZbCAQC/P3vfycSiTB9+nSqqqq4+OKL+de//gVE2+dWVVUxY8YMbrnlFqZMmQLAo48+yiWXXMKFF17I3Llz2bNnD3PmzGHq1KlUVlbGW+hu376dyZMnc/XVVzNlyhQWLlzI2rVrOfPMMznttNNYv359Wq5BKgl9A3CamZWbWSGwAFideIJzrtw5V+acKwNWAT9wzv0+LRGKZMDelsZ4VR4IBAiFQv2eV1wRJBAIANENDVStZ6+//e1v/PCHP2TLli2ccMIJPPHEE/FjH3zwAU899RRbtmyhubmZ22+/nZkzZ3LRRRdx7733EolE+MIXvsBll13GPffcQ3NzM5WVlfz0pz8F4Morr2T58uW8/vrrjB49us/nvv766zz22GO8+OKLFBUV8dRTT/HGG2/w0ksv8aMf/Sje22Xbtm3ceOONNDc389ZbbxEOh3n11Ve57777uOuuu9JyDYYcQ3fOHTCz64jOXhkNPOKc22Jmi3qPL09LJCIjLFaVD6YkUEtj/TIA3y5G8Yvy8vL4D9+vfe1rbN++PX5szJgxFBUVcfXVV3P++edzwQUXHPb63bt38+GHHzJ79mwALr/8ci655BI+/PBDurq6mDlzJgChUKhP24BzzjmHz33uc0C04+Jtt93GK6+8wqhRo3j33Xf55z//GY+vsrISgDPOOIM5c+ZgZlRWVvaJ9WiktFLUObcGWJP0XL+J3Dl3xdGHJSJyZI499tj416NHj+6zM9AxxxzD+vXreeGFF1ixYgUPPfQQL774YkrvO1S/q1hLXoDHH3+cXbt2sXHjRgoKCigrK6O7u/uw+EaNGhV/PGrUqLSNv2vpv+SNhoYGdoaj1XZPRxtMGnvE7xEbT4/pKqqkJFCbrhAlQ/bs2cO+ffs477zzmD59Ol/84hcB+rTAPf744znxxBP54x//yDe/+U1+85vfMHv2bE488URKSkpYt24d06dPZ8WKFQN+zu7duyktLaWgoICXXnqJd955Z0T+fjFK6JI3wuEwPR1tFJaWU1haPuC4+UCSz49EInSP6VRCT5KNnSa7urqYN28e3d3dOOf4+c9/DsCCBQu45pprePDBB1m1ahWPPfYYixYtYt++fUyaNIlf//rXAPzqV7/immuuobi4mGAwyPHHH9/v5yxcuJALL7yQ6upqAoEAkydPHrG/I6h9ruSRYDDIutZOTg7VA4MnnlSaeCW/X6JsTGqZkmvtc4djz549fPaznwWgvr6e999/nwceeCDjn3uk7XNVoYuIDOGZZ57h7rvv5sCBA5x66qk8+uijXofULyV0kSHk68YVcsj8+fOZP3++12EMSQldfCt5SX8kEoExEwd+wTD0dLSxMxxtQFpcEczb8XTn3IhthJwvhjMcrn7o4lvJS/oDgQDFFcG0vX8oFKKwtByIJva9LY1pe+9cUlRURGdn57ASkPTPOUdnZydFRUVH9DpV6OJryYuH0tk4qq6ujrtaxwPEq/R8NGHCBNrb29m1a5fXofhKUVEREyZMOKLXKKGLyFEpKCigvLzc6zAEJXTJIbo5KTI4jaGLiPiEErqIiE9oyEV8JXGqYqzP+UhJnMKoXY7EC0ro4gux8fWd4WUUfbSDQCAwaJ/zdEucDhnb5UgJXUaaErr4Tip9ztOtJFAbX1SUz1MYxVsaQxcR8QkldBERn1BCFxHxCY2hi++sa+2M3yQ92p7nIrlEFbqIiE8ooYuI+ISGXCTnJW/+HGtp66XkzaRDoZDmpUvGqUKXnBfb/BmgsLQ8rT3Ph6O4IthnhWokEumz0YZIpqhCF8+l4+ZkYWl5v5s1p1sq/dRLArU01i+LP06s1EUySRW6iIhPKKGLiPiEErqIiE8ooYuI+IRuikrWSseGzuncFFok26lCFxHxCVXoknO6Is8SDN4bfxyJRGDMRO8CSkHiQiMtMpJMUUKXnLO3pZFI765EEN3QYlNRpbdBDSJx16RIJAKghC4ZoYQuOSl5V6JsHiuvq6uLJ3AtMpJMUkIXyYDkHzBqzysjQTdFRUR8QgldRMQnlNBFRHxCCV1ExCeU0EVEfCKlhG5mtWa21cy2mdnifo7PM7NmM4uYWZOZfSP9oYqIyGCGnLZoZqOBh4FzgHZgg5mtds61JJz2ArDaOefMrAr4T2ByJgIWEZH+pTIPvQbY5pxrBTCzFcA8IJ7QnXN7Es4vBlw6gxTpijzL3pZGILpvKJPGehvQUVAbAMmUVIZcxgM7Eh639z7Xh5ldbGZvAc8AV/X3RmZW1zsk07Rr167hxCt5am9LY599QxOX0+eSUCgUb1mgvUYl3VKp0K2f5w6rwJ1zTwFPmdks4E7g7H7OaQAaAKqrq1XFyxFJ3De0ri43V16qDYBkUioVejuQ2MpuAvDeQCc7514BvmBmJx1lbCIicgRSSegbgNPMrNzMCoEFwOrEE8zsi2ZmvV9PBQqBznQHKyIiAxtyyMU5d8DMrgOeA0YDjzjntpjZot7jy4H/BVxmZp8A+4H5zjkNqYiIjKCUui0659YAa5KeW57w9T3APekNTcSfYp0Yd7Z2Mj1ptk5il0Z1aJQjpZWiIiI+oYQuIuITSugiIj6hhC4i4hPagk6yUkNDAzvDy+KPezraKCwtjz9O9x6iXu1JmtgGAKCrqJKSQO1h52lLO0mFKnTJSuFwOL7UH6KrRIsrgt4FlAHFFcF4GwCIJvdYvxqR4VCFLlkrcam/H5UEammsP/RbSDAYZF2r1uPJ8KlCFxHxCSV0ERGfUEIXEfEJJXQREZ9QQhcR8QkldBERn1BCFxHxCSV0ERGf0MIiyRpdkWcJBu8FoqsmGTNx8Bf4UE9HGzvDiwFomPRufP9RkVSoQpessbelMZrIgUAg4Lul/kMJhULxfjU9HW2Ew2GPI5JcowpdskogEKCxsRHwrmGWV+rq6rirdTxAvEoXORJK6CIeyrcfWpJZGnIREfEJJXQREZ9QQhcR8QkldBERn9BNURlxuhEokhmq0EVEfEIJXUTEJ5TQRUR8QmPo4qmuyLPxne57Otpg0lhvA8oikUiEYDAIQFdRJSWBWm8DkqynhC4jYqAboXtbGunpaKOwtJzC0nJCodAIR5adiiuCVHZvAqKJvXtMpxK6DEkJXTxXWFrOyaF6AOrqzvc4muxQEqilsX4ZAMFgkHWtnR5HJLlAY+giIj6hhC4i4hNK6CIiPqExdJEsFbuRvFPj55IiVegiIj6hhC4i4hNK6CIiPqGELiLiE7opKlnFr611/fr3kuySUkI3s1rgAWA08EvnXH3S8YXArb0P9wDfd879NZ2Bij8k9m4B4sv+ZXA9HW3sDC+OP26Y9C51dXUeRiTZaMghFzMbDTwMnAtUAJeaWUXSaW3AbOdcFXAn0JDuQMUfYr1bYgpLyymuCHoXUA4orgj2+aHX09FGOBz2MCLJVqlU6DXANudcK4CZrQDmAS2xE5xzf0o4fx0wIZ1Bir8k9m6RoZUEavs05kqs1EUSpXJTdDywI+Fxe+9zA/nfwH/1d8DM6sysycyadu3alXqUIiIypFQSuvXznOv3RLNvEU3ot/Z33DnX4Jyrds5Vjxs3LvUoRURkSKkMubQDExMeTwDeSz7JzKqAXwLnOue0VllEZISlUqFvAE4zs3IzKwQWAKsTTzCzU4Ange85595Of5giIjKUISt059wBM7sOeI7otMVHnHNbzGxR7/HlwBJgLPAfZgZwwDlXnbmwRUQkWUrz0J1za4A1Sc8tT/j6auDq9IYmuU6LaURGlpb+i4j4hBK6iIhPqJeLZFzicn8t9U+PSCRCMBgEIBQKqQ2AAKrQZQQkLvfXUv+jV1wRJBAIANHErjYAEqMKXUaElvunT0mglu1EWwF0t6oNgByiCl1ExCeU0EVEfEIJXUTEJ5TQRUR8QgldRMQnlNBFRHxCCV1ExCeU0EVEfEILiyQjGhoa2BleBmi5f6apDYDEqEKXjAiHw1ruPwLUBkASqUKXjNFy/8wrCdTSWB/9TShWpUv+UoUuIuITSugiIj6hhC4i4hMaQ5cjlrhX6Pb68z2MREQSqUIXEfEJJXQREZ9QQhcR8QmNoUtaNDQ09FnUEolEYMxE7wLKU4mrRkErR/ONKnRJi3A4HE3ivQKBgFaHjrBQKBRfNQpaOZqPVKFL2gQCARobG+OPE2fDSObV1dX1qca1cjT/qEIXEfEJJXQREZ9QQhcR8QkldBERn1BCFxHxCSV0ERGfUEKXYeuKPEswGCQYDPaZgy7ZI7bQKBgM0tDQ4HU4kmGahy7DtrelkVdi+4WOmcimokrNPffAQN0vQ6FQ/OvYD1ytGvU3JXQ5KtpmLnslLjTSIqP8oCEXERGfUEIXEfEJJXQREZ/QGLqIjyTflB5si0BtJeg/KVXoZlZrZlvNbJuZLe7n+GQze93MPjazm9MfpoiIDGXICt3MRgMPA+cA7cAGM1vtnGtJOO0D4Abg3zIRpIiIDC2VIZcaYJtzrhXAzFYA84B4QnfOdQAdZqbf23yuoaGBneFlAPTE5qBLTkjezairqJKSQK13AUnapTLkMh7YkfC4vfc5yUPhcJiejjYgOgdduxLlhv52M9rb0uhZPJIZqVTo1s9zbjgfZmZ1QB3AKaecMpy3EA8k3jzb2dqpxUQ5qL/djNa1dnoYkWRCKhV6O5C42+8E4L3hfJhzrsE5V+2cqx43btxw3kJERAaQSkLfAJxmZuVmVggsAFZnNiwRETlSQw65OOcOmNl1wHPAaOAR59wWM1vUe3y5mZ0MNAFjgINmdhNQ4Zz7KHOhi4hIopQWFjnn1gBrkp5bnvD1TqJDMSKSpZLvhYj/aKWoDCpxmiJoqqKf9HS0sTMcXSfYMOldtdb1AfVykUElTlMETVX0i+KKYPwHc09HG+Fw2OOIJB1UocuQNE3Rf0oCtfFFRbEqXXKfKnQREZ9QhS79it1A082z3KYtAfOLKnQREZ9QQhcR8QkNuchh1FEx/yR2YgyFQprCmKNUocth1FExvxRXBOOdGCORiKYw5jBV6NIvTVXMHyWBWrYTncLY3bqYda2dA95M1VZ12U0VuoiITyihi4j4hBK6iIhPaAxdgOjMltjNsEgkAmMmDv4C8a3Epl0QvWmqvUdzgxJ6HuvTTjW8jKKPdhAIBAgEAmwqqvQwMvFK8oym2GwnJfTcoIQucYFAgMbGRkBLxvNVYtMuUOOuXKMxdBERn1BCFxHxCQ255Kn+diJi0lgPI5JspZ2NcocSep6KLe+P9WkpLC0nFAp5HJVkm8SbpD0dbdywdBl3tY4HtGo0Gymh57Hk5f11dfoHKn1pZ6PcojF0ERGfUEIXEfEJDbnkEa0GlaOlG6TZTQk9j4TDYSKRyICrQbWYSAaTfIM0HA4PK6Enfp/pxmp6KaHnGa0GleHSDdLsp4Tuc8nDLLGdaUSOVuK2daCt67KBbor6XGyYBaLVueaaSzokblsH2rouW6hC94nBxiW7x0xk+/RbAM01l/RI3LYOolvXifdUoYuI+IQqdJ/pijxLMHhv/LGmJ8pISRxT13i6N1Sh+8zelsb4mDlEx82TNy0QSbfEMXWNp3tHFboPJHZO7OloY/r06vjURND0RMm8kkAtjfXR78FgMKhq3SNK6Dkslqh3hpfFOycWlpazqahSSVxGXOx7rquoktgEmNhviwMldC0ySi8l9ByVXJUnd04U8cpg1TpEE772KM0MJfQcldjPvLC0XOPkkpWS1z1EIhG6x3QqoWeIEnqOSFzxCdF/GKrKJdvV1dX1GW4JBoO8sq4p3jqguCKo5J5GSuhZLDGJv/zyywDMnj0boN/mWiLZLhQKsa61E4CPd2zm4x2b2dvSCKh7YzqYc86TD66urnZNTU2efHY2a2ho4Ial0fHHj3dsBg4l8U0ae5QckXiDc6Ab9F2RZ+PJPPa9fuzEKfHjD95+fcYSfC7fjDWzjc656v6OqUL3WPJQSqwSP3biFI6dOKXPN7VmroifJHZvTEzuEE3w1157bfzfhqY+pialhG5mtcADwGjgl865+qTj1nv8PGAfcIVz7o00x5qzkpN2ouShlNmzZ/epxNV7RfJBYnKHaIKv7N4ERP+NvPzyywP+G1KyP2TIhG5mo4GHgXOAdmCDma12zrUknHYucFrvn68D/7f3v74zWHIeSHLSTjR79uzDviFViUu+S5z6OFRBNFiyH4hfp06mUqHXANucc60AZrYCmAckJvR5wP9z0QH5dWZ2gpl93jn3/kBvunXr1j5zU3PFYMl5IP0lbRFJTfJMmUTDL7Bejg/xBNfdO+j5uWTIm6Jm9h2g1jl3de/j7wFfd85dl3DO00C9c+7V3scvALc655qS3qsOiP2f+RKwNV1/kaNwEvA/XgeRJXQtDtG1OETX4pBsuBanOufG9XcglQrd+nku+adAKufgnGsAGlL4zBFjZk0D3THON7oWh+haHKJrcUi2X4tUui22A4n9VycA7w3jHBERyaBUEvoG4DQzKzezQmABsDrpnNXAZRY1Hdg92Pi5iIik35BDLs65A2Z2HfAc0WmLjzjntpjZot7jy4E1RKcsbiM6bfHKzIWcdlk1BOQxXYtDdC0O0bU4JKuvhWcrRUVEJL20Y5GIiE8ooYuI+IQSei8zu9nMnJmd5HUsXjGze83sLTNrNrOnzOwEr2MaaWZWa2ZbzWybmS32Oh6vmNlEM3vJzN40sy1mdqPXMXnNzEab2V96191kJSV0ot+8RFsb/MPrWDz238AU51wV8DbwY4/jGVEJbS7OBSqAS82swtuoPHMA+JFz7svAdOCHeXwtYm4E3vQ6iMEooUf9HPg/9LMYKp845553zh3ofbiO6HqCfBJvc+Gc6wFibS7yjnPu/ViDPedcF9FENt7bqLxjZhOA84Ffeh3LYPI+oZvZRcC7zrm/eh1LlrkK+C+vgxhh44EdCY/byeMkFmNmZcBXgT97HIqX7ida9B30OI5B5UU/dDNbC5zcz6F/B24D5o5sRN4Z7Fo45/7Qe86/E/2V+/GRjC0LpNTCIp+Y2WeBJ4CbnHMfeR2PF8zsAqDDObfRzIIehzOovEjozrmz+3vezCqBcuCv0ZbuTADeMLMa59zOEQxxxAx0LWLM7HLgAmCOy79FCmphkcDMCogm88edc096HY+HzgQuMrPzgCJgjJn91jn3XY/jOowWFiUws+1AtXPO625qnujdyORnwGzn3C6v4xlpZnYM0ZvBc4B3iba9CDnntngamAd6N615DPjAOXeTx+Fkjd4K/Wbn3AUeh9KvvB9Dlz4eAkqA/zaziJkt9zqgkdR7QzjW5uJN4D/zMZn3OhP4HnBW7/dCpLdClSymCl1ExCdUoYuI+IQSuoiITyihi4j4hBK6iIhPKKGLiPiEErqIiE8ooYuI+MT/B7rbxyCzFK+SAAAAAElFTkSuQmCC\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"def sample_Zn(alpha,b,n):\n",
|
|
" assert alpha > 2\n",
|
|
" assert b > 0\n",
|
|
" assert n >= 1 and type(n) == int\n",
|
|
" \n",
|
|
" E_X = alpha*b/(alpha - 1)\n",
|
|
" Var_X = alpha*b**2/( (alpha - 1)**2*(alpha - 2) )\n",
|
|
" \n",
|
|
" inv_pareto_samples = [inversion_sample(lambda p: f_inv_pareto(alpha, b, p)) for _ in range(n)]\n",
|
|
" return np.sqrt(n/Var_X)*(np.mean(inv_pareto_samples) - E_X)\n",
|
|
"\n",
|
|
"# Plotting\n",
|
|
"alpha = 4\n",
|
|
"b = 1\n",
|
|
"n = 1000\n",
|
|
"pdf = gaussian\n",
|
|
"samples = [inversion_sample(lambda p: sample_Zn(alpha, b, n)) for _ in range(1000)]\n",
|
|
"compare_plot(samples, pdf, -5, 5, 100)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 10,
|
|
"id": "b5360d77",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "e50b33644ddd6bce391b36cefcc2e308",
|
|
"grade": true,
|
|
"grade_id": "cell-5d16b014bef9d86f",
|
|
"locked": true,
|
|
"points": 10,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"assert_almost_equal(np.mean([sample_Zn(3.5,2.1,100) for _ in range(100)]),0,delta=0.3)\n",
|
|
"assert_almost_equal(np.std([sample_Zn(3.5,2.1,100) for _ in range(100)]),1,delta=0.3)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "6192f05d",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "08ece68d59de21d798d9a955f59be690",
|
|
"grade": false,
|
|
"grade_id": "cell-3e7a23657e9b8374",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"__(b)__ Now take $\\alpha = 3/2$ and $b=1$. \n",
|
|
"With some work (which you do not have to do) one can show that the characteristic function of $X$ admits the following expansion around $t=0$,\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\varphi_X(t) = 1 + 3 i t - (|t|+i t)\\,\\sqrt{2\\pi|t|} + O(t^{2}).\n",
|
|
"$$\n",
|
|
"\n",
|
|
"Based on this, prove the **generalized CLT** for this particular distribution $X$ which states that $Z_n = c\\, n^{1/3} (\\bar{X}_n - \\mathbb{E}[X])$ in the limit $n\\rightarrow\\infty$ converges in distribution, with a to-be-determined choice of overall constant $c$, to a limiting random variable $\\mathcal{S}$ with characteristic function \n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\varphi_{\\mathcal{S}}(t) = \\exp\\big(-(|t|+it)\\sqrt{|t|}\\big).\n",
|
|
"$$\n",
|
|
"\n",
|
|
"__(15 pts)__"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "9735cd88",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "dfd8683eea5663baa81f138a2809722b",
|
|
"grade": true,
|
|
"grade_id": "cell-b25551eca32c4807",
|
|
"locked": false,
|
|
"points": 15,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"\\begin{align}\n",
|
|
"\\phi_{Z_n}(t) &= \\mathbb{E}\\left[ e^{itZ_n} \\right]\n",
|
|
" \\\\ &= \\mathbb{E}\\left[ e^{itcn^{1/3}(\\bar{X_n} - \\mathbb{E}[X])} \\right]\n",
|
|
" \\\\ &= \\mathbb{E}\\left[ e^{itcn^{1/3}(\\frac{1}{n}\\sum_{i=1}^n X_i - \\mathbb{E}[X])} \\right]\n",
|
|
" \\\\ &= \\mathbb{E}\\left[ \\left( \\prod_{i=1}^n e^{itcn^{-2/3}X_i} \\right) e^{itcn^{1/3}\\mathbb{E}[X])} \\right]\n",
|
|
" \\\\ &= \\left( \\prod_{i=1}^n \\mathbb{E}\\left[ e^{itcn^{-2/3}X_i} \\right] \\right)\\mathbb{E}\\left[ e^{itcn^{1/3}\\mathbb{E}[X])} \\right]\n",
|
|
" \\\\ &= \\left( \\prod_{i=1}^n \\phi_X(cn^{-2/3}t) \\right)\\mathbb{E}\\left[ e^{itcn^{1/3}\\mathbb{E}[X])} \\right]\n",
|
|
" \\\\ &= \\left( \\phi_X(cn^{-2/3}t) \\right)^n \\mathbb{E}\\left[ e^{itcn^{1/3}\\mathbb{E}[X])} \\right]\n",
|
|
"\\end{align}\n",
|
|
"where we used the identity for products of indepedent expectation values <https://hef.ru.nl/~tbudd/mct/lectures/probability_random_variables.html#equation-product-expectation>, and the definition of $\\phi_X(t) := \\mathbb{E}\\left[ e^{itX} \\right]$.\n",
|
|
"\n",
|
|
"Next, we will use the Taylor expansion around $t = 0$ as is given above, and, for the latter exponential, $\\mathbb{E}(X) = 3$ for $\\alpha = 3/2, b = 1$ as given.\n",
|
|
"\n",
|
|
"\\begin{align}\n",
|
|
"\\phi_{Z_n}(t) &= \\left( \\phi_X(cn^{-2/3}t) \\right)^n \\mathbb{E}\\left[ e^{itcn^{1/3}\\mathbb{E}[X])} \\right]\n",
|
|
" \\\\ &= \\left( 1 + 3 i cn^{-2/3}t - (|cn^{-2/3}t|+i cn^{-2/3}t)\\,\\sqrt{2\\pi|cn^{-2/3}t|} + \\mathcal{O}(t^2) \\right)^n e^{3itcn^{1/3}}\n",
|
|
" \\\\ &= \\left( 1 + \\frac{1}{n} \\left[ 3 i cn^{1/3}t - (|ct|+i ct)\\,\\sqrt{2\\pi|ct|} \\right] + \\mathcal{O}(t^2) \\right)^n e^{3itcn^{1/3}}\n",
|
|
"\\end{align}\n",
|
|
"\n",
|
|
"Taking the limit $n \\to \\infty$, the first set of parentheses can be rewritten in terms of an exponential using the identity $\\lim_{n\\to\\infty} (1 + \\frac{a}{n})^n = e^{a}$, we find a way to our desired expression.\n",
|
|
"\n",
|
|
"\\begin{align}\n",
|
|
"\\lim_{n\\to\\infty} \\phi_{Z_n}(t) &= \\lim_{n\\to\\infty} \\left( 1 + \\frac{1}{n} \\left[ 3 i cn^{1/3}t - (|ct|+i ct)\\,\\sqrt{2\\pi|ct|} \\right] + \\mathcal{O}(t^2) \\right)^n e^{-3itcn^{1/3}}\n",
|
|
" \\\\ &= \\lim_{n\\to\\infty} \\exp{({3 i cn^{1/3}t - (|ct|+i ct)\\,\\sqrt{2\\pi|ct|}})} \\exp{(e^{-3itcn^{1/3}})}\n",
|
|
" \\\\ &= \\exp{({-(|ct|+i ct)\\,\\sqrt{2\\pi|ct|}})}\n",
|
|
"\\end{align}\n",
|
|
"\n",
|
|
"This matches $\\phi_S(t) = \\exp\\big(-(|t|+it)\\sqrt{|t|}\\big)$ for $\\sqrt{2\\pi} c^{3/2} = 1 \\implies c = (2\\pi)^{\\frac{-1}{3}}$."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "5b1d9f54",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "342020128f929d47eabfdf9c075ff20c",
|
|
"grade": false,
|
|
"grade_id": "cell-d1701433c3c77172",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"__(c)__ The random variable $\\mathcal{S}$ has a [stable Lévy distribution](https://en.wikipedia.org/wiki/Stable_distribution) with index $\\alpha = 3/2$ and skewness $\\beta = 1$. Its probability density function $f_{\\mathcal{S}}(x)$ does not admit a simple expression, but can be accessed numerically using SciPy's `scipy.stats.levy_stable.pdf(x,1.5,1.0)`. Verify numerically that the generalized CLT of part (b) holds by comparing an appropriate histogram to this PDF. __(10 pts)__"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 11,
|
|
"id": "b06896e5",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "c6fe081427f342c354ee8a9b3b3331e7",
|
|
"grade": true,
|
|
"grade_id": "cell-e08d054985cfa762",
|
|
"locked": false,
|
|
"points": 10,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD7CAYAAABkO19ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAg0UlEQVR4nO3de3DU5dn/8fcFYlNTUIugHQ6KLdZGErc0IkgfWatSPPJzOg4YW619NNJ6qDPVX6mPdabWqTh12qr1eTA/a/VpTaWjYhmlaq0GqpVCsFsEFEsDlogWGk8BiYhcvz/2wDfLJvkm2c1uvvt5zTDke9rcu8qVO9d939dt7o6IiETXkGI3QERECkuBXkQk4hToRUQiToFeRCTiFOhFRCJOgV5EJOJCBXozm2VmG8xso5nNz3F9tpmtMbOEmTWb2RcD1zab2Uvpa/lsvIiI9Mx6mkdvZkOBV4HTgVZgFXCBu68P3PMJYKe7u5nVAL9192NT1zYDte7+78K8BRER6c4BIe6ZAmx09xYAM3sQmA1kAr277wjcXwn0axXWYYcd5kcddVR/XkJEpKysXr363+4+Kte1MIF+DLAlcNwKnJh9k5mdB9wCjAbOClxy4Ckzc+Bud2/I9U3MrB6oBxg/fjzNzcryiIiEZWavdXUtTI7ecpzbr8fu7otT6Zr/A/wwcGm6u08GzgCuMLOTc30Td29w91p3rx01KucPJRER6YMwgb4VGBc4Hgts7epmd18OfNrMDksdb039vQ1YTDIVJCIiAyRMoF8FTDSzCWZ2IDAXWBK8wcw+Y2aW+noycCDQZmaVZjY8db4SmAmszecbEBGR7vWYo3f3PWZ2JfAkMBS4193Xmdm81PWFwFeAi8zsQ2AXMCc1A+dwYHHqZ8ABQKO7P1Gg9yIiJebDDz+ktbWVjo6OYjclMioqKhg7dizDhg0L/UyP0yuLoba21jUYKzL4bdq0ieHDhzNy5EhSHT7pB3enra2N9vZ2JkyY0Omama1299pcz2llrIgUTEdHh4J8HpkZI0eO7PVvSAr0IlJQCvL51ZfPU4FeRCTiFOhFRPIokUiwdOnSXj8Xj8cLtlA0zMpYkaI4av7jma83LzirmztFSkcikaC5uZkzzzyz2E3JUI9eRCLt17/+NVOmTCEWi3H55Zfzl7/8hZqaGjo6Oti5cyfHHXcca9eupampiZNPPpnzzjuPqqoq5s2bx969ewF46qmnmDZtGpMnT+b8889nx45kea9Vq1Zx0kkncfzxxzNlyhTeffddbrzxRhYtWkQsFmPRokXs3LmTb3zjG5xwwgl8/vOf53e/+x0Au3btYu7cudTU1DBnzhx27dpVsM9APXoRGRDXXHMNiUQir68Zi8X42c9+1uX1l19+mUWLFvH8888zbNgwvvWtb7FhwwbOPfdcbrjhBnbt2sVXv/pVJk2aRFNTEytXrmT9+vUceeSRzJo1i0ceeYR4PM7NN9/M008/TWVlJbfeeis/+clPmD9/PnPmzGHRokWccMIJvPfeexx00EHcdNNNNDc38/Of/xyA66+/ni996Uvce++9vPPOO0yZMoXTTjuNu+++m4MOOog1a9awZs0aJk+enNfPJkiBXgoimHYBpV6kOP74xz+yevVqTjjhBCDZix49ejQ33ngjJ5xwAhUVFdxxxx2Z+6dMmcLRRx8NwAUXXMBzzz1HRUUF69evZ/r06QDs3r2badOmsWHDBj71qU9lXnvEiBE52/DUU0+xZMkSbrvtNiA55fSf//wny5cv5+qrrwagpqaGmpqawnwIKNCLyADpruddKO7OxRdfzC233NLp/JtvvsmOHTv48MMP6ejooLKyEth/6qKZ4e6cfvrp/OY3v+l0bc2aNaGmOro7Dz/8MJ/97Gf3uzZQU0+VoxeRyDr11FN56KGH2LZtGwBvvfUWr732GvX19fzwhz/kwgsv5Lvf/W7m/pUrV7Jp0yb27t3LokWL+OIXv8jUqVN5/vnn2bhxIwDvv/8+r776Ksceeyxbt25l1apVALS3t7Nnzx6GDx9Oe3t75jW//OUvc+edd5KuQvDXv/4VgJNPPpkHHngAgLVr17JmzZqCfQ7q0YtIZFVVVXHzzTczc+ZM9u7dy7Bhw5g9ezYHHHAAdXV1fPTRR5x00kk888wzDBkyhGnTpjF//nxeeumlzMDskCFDuO+++7jgggv44IMPALj55ps55phjWLRoEVdddRW7du3i4x//OE8//TSnnHIKCxYsIBaL8b3vfY/vf//7XHPNNdTU1ODuHHXUUTz22GN885vf5JJLLqGmpoZYLMaUKYUr7KtaN1IQ+cjRZ79Gf19PBt7LL7/M5z73uWI3I5SmpiZuu+02HnvssWI3pUe5PlfVuhERKWNK3YiIkFyZGo/Hi92MglCPXkQk4hToRUQiToFeRCTiFOhFRCJOg7EiMmC6mzLbF2Gm2W7evJmzzz6btWvXdjp/4403cvLJJ3PaaaflfO7RRx/lmGOOoaqqKi9tLSb16EWkLN10001dBnlIBvr169fn5Xvt2bMnL6/TVwr0IhJ5H330EZdddhnHHXccM2fOZNeuXXz961/noYceAmD+/PlUVVVRU1PDtddey5///GeWLFnCddddRywW4x//+AeJRIKpU6dSU1PDeeedx9tvvw0kSxXX1NQwbdo0rrvuOiZNmgTAfffdx/nnn88555zDzJkz2bFjB6eeeiqTJ0+muro6U6548+bNHHvssVx66aVMmjSJCy+8kKeffprp06czceJEVq5c2e/3r0AvIpH397//nSuuuIJ169ZxyCGH8PDDD2euvfXWWyxevJh169axZs0abrjhBk466STOPfdcfvzjH5NIJPj0pz/NRRddxK233sqaNWuorq7mBz/4AQCXXHIJCxcu5IUXXmDo0KGdvu8LL7zA/fffzzPPPENFRQWLFy/mxRdf5Nlnn+U73/lOpv7Nxo0b+fa3v82aNWt45ZVXaGxs5LnnnuO2227jRz/6Ub/ff6hAb2azzGyDmW00s/k5rs82szVmljCzZjP7YthnpTwcNf/xzJ9SfD2JtgkTJhCLxQD4whe+wObNmzPXRowYQUVFBZdeeimPPPIIBx100H7Pv/vuu7zzzjvMmDEDgIsvvpjly5fzzjvv0N7ezkknnQRAXV1dp+dOP/10PvnJTwLJKpbXX389NTU1nHbaabz++uv861//yrSvurqaIUOGcNxxx3HqqadiZlRXV3dqa1/1GOjNbChwF3AGUAVcYGbZoxN/BI539xjwDeCeXjwrIlJQH/vYxzJfDx06tFPO/IADDmDlypV85Stf4dFHH2XWrFmhX7enWmHp8scADzzwANu3b2f16tUkEgkOP/xwOjo69mvfkCFDMsdDhgzJS34/zKybKcBGd28BMLMHgdlAZpTC3XcE7q8EPOyzUn60F6yUkh07dvD+++9z5plnMnXqVD7zmc8AdCo3fPDBB3PooYfypz/9if/4j//gV7/6FTNmzODQQw9l+PDhrFixgqlTp/Lggw92+X3effddRo8ezbBhw3j22Wd57bXXBuT9QbhAPwbYEjhuBU7MvsnMzgNuAUYD6X+9oZ5NPV8P1AOMHz8+RLNEZLApxR/s7e3tzJ49m46ODtydn/70pwDMnTuXyy67jDvuuIOHHnqI+++/n3nz5vH+++9z9NFH88tf/hKAX/ziF1x22WVUVlYSj8c5+OCDc36fCy+8kHPOOYfa2lpisRjHHnvsgL3HHssUm9n5wJfd/dLU8deAKe5+VRf3nwzc6O6n9fbZNJUpHvzC5s67+4efj9eQ4hpMZYr7aseOHXziE58AYMGCBbzxxhvcfvvtBf2evS1THKZH3wqMCxyPBbZ2dbO7LzezT5vZYb19VqQvtD+tFNPjjz/OLbfcwp49ezjyyCO57777it2k/YQJ9KuAiWY2AXgdmAt0Glo2s88A/3B3N7PJwIFAG/BOT8+KiAxmc+bMYc6cOcVuRrd6DPTuvsfMrgSeBIYC97r7OjObl7q+EPgKcJGZfQjsAuZ4MieU89kCvRcpMk11lFzcfcA2wS4HfdkVMFStG3dfCizNOrcw8PWtwK1hn5XoUHCX7lRUVNDW1sbIkSMV7PPA3Wlra6OioqJXz6momRSV8uvRNnbsWFpbW9m+fXuxmxIZFRUVjB07tlfPKNCLSMEMGzaMCRMmFLsZZU+1bkREIk49eikpyvmL5J969CIiEadALyIScUrdSMloTzzBzvVNmePKqjjDY+ErCYpIbgr0UjTZgf2DLck9PT82bhIfbFnLB1vWZq4r6Iv0nQK9DKhgcA8G9vTf6YAevG/3tk0ACvQifaRALwNq5/omdm/bxIGjJ3QK7NmGx2Zlzr/ZOJ/d2zbxZuO+DcrUwxcJT4FeBtyBoydwRN2C0PdXVsU7HauHL9I7CvRScNlpmANH926lZLB3D/v38NW7F+mepldKwaXTNZDszWf30Hursiqe+WGxe9umTgO6IrI/9ehlQPQ2XdOd7Py9iHRPgV4iR5uPi3SmQC95lz0/vi95+d7QjByR7ilHL3kXzMlDfvLyXQnm60E5e5Fc1KOXgshnTr47uWbkiEhn6tGLiEScevQSOcGcfcPRr1NfX1/kFokUlwK95EV/F0XlS3AsYPe2TTQ2NirQS9lToJe8CNawKeTga0+y59ivaGnrctcqTb2UchEq0JvZLOB2YChwj7svyLp+IfDd1OEO4Jvu/rfUtc1AO/ARsMfda/PTdCk1AzUAKyK902OgN7OhwF3A6UArsMrMlrj7+sBtm4AZ7v62mZ0BNAAnBq6f4u7/zmO7RUQkpDA9+inARndvATCzB4HZQCbQu/ufA/evAMbms5EifaXiZyLhpleOAbYEjltT57ryn8DvA8cOPGVmq82sy1ExM6s3s2Yza96+fXuIZol0T8XPRJLC9OgtxznPeaPZKSQD/RcDp6e7+1YzGw38wcxecffl+72gewPJlA+1tbU5X1+KLziwWSozbbqi4mciSWF69K3AuMDxWGBr9k1mVgPcA8x297b0eXffmvp7G7CYZCpIIiDf5YdFpDDC9OhXARPNbALwOjAXqAveYGbjgUeAr7n7q4HzlcAQd29PfT0TuClfjZfi00wbkdLXY6B39z1mdiXwJMnplfe6+zozm5e6vhC4ERgJ/LeZwb5plIcDi1PnDgAa3f2JgrwTkR5kV7nUqlkpF6Hm0bv7UmBp1rmFga8vBS7N8VwLcHw/2yjSb7n2ndWqWSkXWhkrZUFVLqWcqXqliEjEKdCLiEScUjcSWkNDA2823pk5LsW58yKyPwV6Ca2xsbFTcB/sc+cTiQTxeByAlyqqMzl8VbWUqFGgl5yCK2CDgS8q8+Yrq+LsXN/Eipa21A+vNtXBkchSoJeypPIIUk40GCsiEnEK9CIiEadALyISccrRi9C5Do5q4EjUqEcv3WpPPEE8Hicej5NIJIrdnILI3qCksbGxyC0SyS/16KVbO9c3kXhvC7FYjFgsxksV1cVuUt5pBo5EnQK99KhjxDg2T70OgOFFbstAWNHS1uU6ApHBSKkbEZGIU6AXEYk4pW5EsmgnKoka9ehFAoIzcECzcCQa1KMXCdBOVBJF6tGL9CA9Cyc4E0dkMFGPXvYT3GBEm4uIDH7q0ct+0huMwODfXEREQvbozWwWcDswFLjH3RdkXb8Q+G7qcAfwTXf/W5hnpTRFZYMREQkR6M1sKHAXcDrQCqwysyXuvj5w2yZghru/bWZnAA3AiSGfFSlpKngmg12Y1M0UYKO7t7j7buBBYHbwBnf/s7u/nTpcAYwN+6xIKVPBM4mCMKmbMcCWwHErcGI39/8n8PvePmtm9UA9wPjx40M0S6TwVPBMoiBMj95ynPOcN5qdQjLQp/P1oZ919wZ3r3X32lGjRoVoloiIhBGmR98KjAscjwW2Zt9kZjXAPcAZ7t7Wm2dFRKRwwvToVwETzWyCmR0IzAWWBG8ws/HAI8DX3P3V3jwrIiKF1WOgd/c9wJXAk8DLwG/dfZ2ZzTOzeanbbgRGAv9tZgkza+7u2QK8D+mHhoYGKsZXZ/4sX9Fc7CaVrEQikdlxq6GhodjNEQnF3HOmzIuqtrbWm5sVbAZKPB5n+YrmTitgK6vinWq+SHJbxeqOl4BkwI/FYjQ1NRW3USIpZrba3WtzXVMJBAG0QCqM4bFZbCb5w6+jRTNwZPBQCQQRkYhToBcRiTgFehGRiFOOXqSP0jNw0urq6lQHR0qSAr1IH1RWxTMzcCAZ9AEFeilJCvQifTA8NoumBXdmjoM9e5FSoxy9iEjEqUdfphoaGjIldxOJBIwY1/0DIjJoqUdfphobGzN55Vgspu0CRSJMPfoyFlzCf9T8x4vbGBEpGAX6MraipU0BPo+C0y011VJKiQK9SB7U1dVlvtZUSyk1CvQieVBfX58J7JpqKaVGg7EiIhGnQC9SAOnxD42BSClQoBcRiTjl6EUKYPe2TbzZmNycpOHo1zUwK0WlQF9GtBo2v7LTMpsXnAUkZ+CsaGkDkgG/sbFRgV6KSoG+jDQ2Nu7bG3bEOK2GLZD6+np+1DIGINOrFykmBfoyo71hRcqPBmNFRCIuVKA3s1lmtsHMNprZfr+LmtmxZvaCmX1gZtdmXdtsZi+ZWcLMmvPVcBERCafHQG9mQ4G7gDOAKuACM6vKuu0t4Grgti5e5hR3j7l7bX8aKzIYLV/RTMX46syfhoaGYjdJykyYHP0UYKO7twCY2YPAbGB9+gZ33wZsM7OzCtJKkUEg1+Ko7AFvzcKRYggT6McAWwLHrcCJvfgeDjxlZg7c7e45uzNmVg/UA4wfP74XLy9SuobHZjE8NitzrFk4UgxhcvSW45z34ntMd/fJJFM/V5jZyblucvcGd69199pRo0b14uVFBheVR5CBFibQtwLBlTVjga1hv4G7b039vQ1YTDIVJCIiAyRMoF8FTDSzCWZ2IDAXWBLmxc2s0syGp78GZgJr+9pYkShIl0d4s3G+BmZlQPSYo3f3PWZ2JfAkMBS4193Xmdm81PWFZnYE0AyMAPaa2TUkZ+gcBiw2s/T3anT3JwryTmQ/wZIHoLIHpSA4OKuBWRko5t6bdPvAqK2t9eZmTbnvr3g8TiKRIBaLZc69VFHdaXBQiufNxvlMPXpkZt9ekf4ws9VdTWFXCYSIC24ADtoEXKQcqQSCiEjEqUcfcempfCJSvtSjFymiRCJBPB4nHo9rBo4UjHr0IkVSWRWnuuMlIDUjCjQDRwpCgV6kSIbHZtG04E4gOUNKpFCUuhERiTgFehGRiFPqJmK0AfjglR6YTaurq1POXvJCPfqIaWxszAzsxWIxbQA+SNTV1XVawZxIJDqVrxDpD/XoIyi4GlZz6AeH+vr6Tr13Dc5KPqlHLyIScQr0IiIRp0AvUqK0albyRTl6kSIKjqFsXnBW5uu6ujpWtLSxoqWN3ds2AVo1K32nQC9Sgurr6/lRyxhAG4pL/yl1IyIScQr0EdDQ0JDJ5abn0IuIpCl1EwHpRVKxWIxYLEZdXV2xmyR90N2ah+UrmqkYXw3AHTdcpXy99IoCfURkbxko0aENxaW/FOhFStzw2KzMhu4amJW+UI4+ItJbBqrkgYhkCxXozWyWmW0ws41mtl+XwsyONbMXzOwDM7u2N8+KSO8EF1JpMZWE0WOgN7OhwF3AGUAVcIGZVWXd9hZwNXBbH54VkZAqq+Kqcim9FiZHPwXY6O4tAGb2IDAbWJ++wd23AdvM7KzePiv5p/RNdAW3HwRVuZRwwqRuxgBbAsetqXNhhH7WzOrNrNnMmrdv3x7y5UVEpCdhAr3lOOchXz/0s+7e4O617l47atSokC8vIiI9CZO6aQWC+9GNBbaGfP3+PCtdCG4XCNoysNwFtyDU9oOSS5ge/SpgoplNMLMDgbnAkpCv359npQvB7QJBWwaWs+AWhBqYla702KN39z1mdiXwJDAUuNfd15nZvNT1hWZ2BNAMjAD2mtk1QJW7v5fr2QK9l7KSvRJWA7DlKVjlsqNFs5clt1ArY919KbA069zCwNdvkkzLhHpWRPquux/qSuNILiqBIBIRlVVxqjteAsik9hToBRToRSIjOMde8+slSLVuREQiTj16kYgK5utBOftypkA/SATnzqc3GRHpSvbmM8rZlzcF+kFCu0hJb9TX13cK6tk5++yZO5sXZJepkihRoB9EtIuU9CQYwBW8JU2BXqRMBHP27RXVmV2rJPoU6EUiKti7b6+oJj2sk0gk6BjRpkBfRjS9UqQMDI/NoqmpiaamJg3klyH16EXK0O5tmzptNN5w9OuakRNhCvSDlIqYSV/V1dWxoqUtc7x72yYaGxsV6CNMgb6Eae68FEKw4iXQqWcv0aQcfQkL1p3X3HkppPSMnHg8TkNDQ7GbI3mmHn2J09x5KTRVvYw+BXqRMqeql9Gn1I2IdKI0TvSoRz+IaKaN9Ed3//+kr7VXVNMxoo0VLW3s3rYJUBonChToRSRjeGxWZsXsm43zVeo4IhToS0hwOiVoSqUUV3CQFjRQO5gp0JeQYCli0JRKKa7gIC1ooHYwU6AvMZpOKaVs+YpmKsZXA3DHDVepdz9IhAr0ZjYLuB0YCtzj7guyrlvq+pnA+8DX3f3F1LXNQDvwEbDH3Wvz1noRKajsCpgHjk6WTlDZhMGlx0BvZkOBu4DTgVZglZktcff1gdvOACam/pwI/E/q77RT3P3feWu1iAy47gZqNUhb2sLMo58CbHT3FnffDTwIzM66Zzbwv560AjjEzD6V57aKSImorIpnxpISiUSnSQRSesKkbsYAWwLHrXTurXd1zxjgDcCBp8zMgbvdPecKDDOrB+oBxo8fH6rxIlIc2atpNQ2ztIUJ9JbjnPfinunuvtXMRgN/MLNX3H35fjcnfwA0ANTW1ma/fmSpQqUMdtkzwzQNs/SECfStwLjA8Vhga9h73D399zYzW0wyFbRfoC9XwSmVmk4pg8m+gdoxMPW6zPmOFuXvS02YQL8KmGhmE4DXgblAdjRaAlxpZg+STOu86+5vmFklMMTd21NfzwRuyl/zo0FTKiVKVA2z9PQ4GOvue4ArgSeBl4Hfuvs6M5tnZvNSty0FWoCNwP8DvpU6fzjwnJn9DVgJPO7uT+T5PYhICcnenzY9975ifLWKpBVJqHn07r6UZDAPnlsY+NqBK3I81wIc3882isggFdy28IMta7n88ss7zdBRWmdgaGXsAOupnk12hcHNC84aqKaJ5F1w28L2xBOqnVMkCvQDTPVspFzlqp2jQduBoUBfBL0ZfFUNeomqYAdn2bJlLFu2LPPbroJ+finQi0jehemg1NfXZ4J5Q0MDV998pzY8KRAFehEpumAuXxue5J8C/QDQ6leRpDCTDbI3PFFap/8U6AeAVr+K5JYrxZM9aBvsKCno940C/QDR6leRvgmmdT5ZUZ3p7WcHfVDg70qYMsUiIiUhuOr27rvvZsaMGZlry5Yt4/LLLycejxOPx7UKN0A9+gJRXl6kMILF1DYHfkvOleK5+uZkCqjctz1UoC8Q5eVFBlZ2imfn+iZApRdAgb6glJcXKayu5usHtz3MLr2QndtP1+JJi2LvX4E+T3qqYSMixdHdLJ5s2b3/qPT8LVl4srTU1tZ6c3NzsZvRK+m6HcHg3pf/SVTyQCT/uisOGPw3F+z9L1u2DKDTgG9Qqf0QMLPV7l6b65p69HnU11SNgrtIaQj2/rN7/sEUzwdb1u43tTPtpYrqTNoISqMCrQJ9P/R1Zo0Cu8jA6su/ueDgLsARU/ddy877pyV/C1iWGQgGiK/4cebrYv0WoNRNP2Sna8L+R1SgFxn8cvXU08XZcvlgy1qg61RQUPC3grC/EXSXulGg76VcvfjepmsU6EWiJRiMu/r33Z54olNPf+rRI3Pelx4b+Ni4Sd3eB507l8rR55Hmx4tItjCdt+CUT4CmLnrq3f1WEJSrBERXFOh70NW0SfXiRaQ/uo4JYziibkHmKPsHQvq5YN2fnijQ55C9lBr25dXUixeRgdTdorDN7PsNgVSsykWBPqWr4D5jxoySmy8rItIboQK9mc0CbgeGAve4+4Ks65a6fibwPvB1d38xzLPFkp2SyUdwD7OpgojIQOsx0JvZUOAu4HSgFVhlZkvcfX3gtjOAiak/JwL/A5wY8tl+625Jc1eyUzLquYtIVIXp0U8BNrp7C4CZPQjMBoLBejbwv56cq7nCzA4xs08BR4V4dj8bNmzotF9kT3paqpxLfwJ72IFVDcCKSCkIE+jHAFsCx60ke+093TMm5LMAmFk9kI66O5YtW7YhRNs6WdbNYERX919++eVdXT4M+Hdv2xBR+iz20Wexjz6LfUrhsziyqwthAr3lOJe9yqqre8I8mzzp3gCUzJYwZtbc1eKDcqPPYh99Fvvos9in1D+LMIG+FRgXOB4LbA15z4EhnhURkQIKs2fsKmCimU0wswOBucCSrHuWABdZ0lTgXXd/I+SzIiJSQD326N19j5ldCTxJcorkve6+zszmpa4vBJaSnFq5keT0yku6e7Yg7yT/SiaNVAL0Weyjz2IffRb7lPRnUZJFzUREJH/CpG5ERGQQU6AXEYk4BfoQzOxaM3MzO6zYbSkWM/uxmb1iZmvMbLGZHVLsNg0kM5tlZhvMbKOZzS92e4rFzMaZ2bNm9rKZrTOzbxe7TcVmZkPN7K9m9lix29IVBfoemNk4kiUc/lnsthTZH4BJ7l4DvAp8r8jtGTCBUh5nAFXABWZWVdxWFc0e4Dvu/jlgKnBFGX8Wad8GXi52I7qjQN+znwL/ly4WepULd3/K3fekDleQXBNRLjJlQNx9N5Au5VF23P2NdMFCd28nGeDGdP9UdJnZWOAs4J5it6U7CvTdMLNzgdfd/W/FbkuJ+Qbw+2I3YgB1VeKjrJnZUcDngb8UuSnF9DOSHcG9RW5Ht8q+Hr2ZPQ0ckePSfwHXAzMHtkXF091n4e6/S93zXyR/fX9gINtWZKFLeZQLM/sE8DBwjbu/V+z2FIOZnQ1sc/fVZhYvcnO6VfaB3t1Py3XezKqBCcDfkuX2GQu8aGZT3P3NAWzigOnqs0gzs4uBs4FTvbwWYIQpA1I2zGwYySD/gLs/Uuz2FNF04FwzOxOoAEaY2a/d/atFbtd+tGAqJDPbDNS6e7Er1BVFagOZnwAz3H17sdszkMzsAJID0KcCr5Ms7VE3iFZ5501qk6H7gbfc/ZoiN6dkpHr017r72UVuSk7K0UtYPweGA38ws4SZLSx2gwZKahA6XcrjZeC35RjkU6YDXwO+lPr/IJHq0UoJU49eRCTi1KMXEYk4BXoRkYhToBcRiTgFehGRiFOgFxGJOAV6EZGIU6AXEYm4/w/lz2dwZIAQ2wAAAABJRU5ErkJggg==\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"from scipy.stats import levy_stable\n",
|
|
"\n",
|
|
"def sample_Zn(alpha, beta, c, n):\n",
|
|
" assert n >= 1 and type(n) == int\n",
|
|
" \n",
|
|
" E_X = alpha*b/(alpha - 1)\n",
|
|
" \n",
|
|
" samples = [inversion_sample(lambda p: f_inv_pareto(alpha, beta, p)) for _ in range(n)]\n",
|
|
" return c*n**(1/3)*(np.mean(samples) - E_X)\n",
|
|
"\n",
|
|
"alpha = 3/2\n",
|
|
"beta = 1\n",
|
|
"c = (2*np.pi)**(-1/3)\n",
|
|
"n = 1000\n",
|
|
"pdf = lambda x: levy_stable.pdf(x, alpha, beta)\n",
|
|
"samples = [sample_Zn(alpha, b, c, n) for _ in range(10000)]\n",
|
|
"compare_plot(samples, pdf, -5, 5, 100)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "f49856d8",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "d8c57e5a527eaad8318e7d31dba01694",
|
|
"grade": false,
|
|
"grade_id": "cell-bc80caacda124bf9",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"## Joint probability density functions and sampling the normal distribution \n",
|
|
"__(30 Points)__\n",
|
|
"\n",
|
|
"Let $\\Phi$ be a uniform random variable on $(0,2\\pi)$ and $R$ an independent continuous random variable with probability density function $f_R(r) = r\\,e^{-r^2/2}$ for $r>0$. Set $X = R \\cos \\Phi$ and $Y = R \\sin \\Phi$. This is called the **Box-Muller transform**.\n",
|
|
"\n",
|
|
"__(a)__ Since $\\Phi$ and $R$ are independent, the joint probability density of $\\Phi$ and $R$ is $f_{\\Phi,R}(\\phi,r) = f_\\Phi(\\phi)f_R(r) = \\frac{1}{2\\pi}\\, r\\,e^{-r^2/2}$. Show by change of variables that $X$ and $Y$ are also independent and both distributed as a standard normal distribution $\\mathcal{N}$. __(15 pts)__"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "aa3821de",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "2514e6664aeb4e24a9e881522a8f3a0f",
|
|
"grade": true,
|
|
"grade_id": "cell-4f20e3b730ba0d23",
|
|
"locked": false,
|
|
"points": 15,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"The coordinate transformation $T$ is defined as $T(\\Phi, R) = (R\\cos{\\Phi}, R\\sin{\\Phi}) = (X, Y)$. As $T$ is invertible differentiable, we can write the equality between the joint probability density in both coordinate pairs as follows, using the Jacobian.\n",
|
|
"$$\n",
|
|
"f_{X,Y}(x,y) \\Big|\\frac{\\mathrm{d}x}{\\mathrm{d}\\phi}\\frac{\\mathrm{d}y}{\\mathrm{d}r}-\\frac{\\mathrm{d}y}{\\mathrm{d}\\phi}\\frac{\\mathrm{d}x}{\\mathrm{d}r}\\Big|\n",
|
|
"= f_{X,Y}(T(\\phi,r)) \\Big|\\frac{\\mathrm{d}x}{\\mathrm{d}\\phi}\\frac{\\mathrm{d}y}{\\mathrm{d}r}-\\frac{\\mathrm{d}y}{\\mathrm{d}\\phi}\\frac{\\mathrm{d}x}{\\mathrm{d}r}\\Big|\n",
|
|
"= f_{X,Y}(T(\\phi,r)) \\Big|-r\\sin{\\phi}\\sin{\\phi}-r\\cos{\\phi}\\cos{\\phi}\\Big|\n",
|
|
"= f_{X,Y}(T(\\phi,r)) r\n",
|
|
"= f_{\\Phi,R}(\\phi,r)\n",
|
|
"= \\frac{1}{2\\pi}\\, r\\,e^{-r^2/2}\n",
|
|
"\\\\ \\implies f_{X,Y}(x,y) = \\frac{1}{2\\pi}\\,e^{-r^2/2} = \\frac{1}{\\sqrt{2\\pi}}e^{-x^2/2}\\frac{1}{\\sqrt{2\\pi}}e^{-y^2/2} = f_{X}(x)f_Y(y)\n",
|
|
"$$\n",
|
|
"using that $r^2 = x^2 + y^2$ in the second to last step. We conclude that $X$ and $Y$ are independent, and the factorization shows that they are both distributed as a standard normal distribution $\\mathcal{N}(0,1)$."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "5d064cef",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "1af73334332fe512ef7d0edb5803a58d",
|
|
"grade": false,
|
|
"grade_id": "cell-2f07fdb2a906bb71",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"__(b)__ Write a function to sample a pair of independent normal random variables using the Box-Muller transform. Hint: to sample $R$ you can use the inversion method of the first exercise. Produce a histogram to check the distribution of your normal variables. __(15 pts)__"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "36710a33",
|
|
"metadata": {},
|
|
"source": [
|
|
"For the sampling of $R$, we take its PDF, calculate its CDF, invert it, and use the function `inversion_sample` to pull values for $R$.\n",
|
|
"\n",
|
|
"$$\n",
|
|
"f_R(r) = re^{-r^2/2}\n",
|
|
"\\\\ \\implies F_R(r) = \\int_{0}^r te^{-t^2/2}dt = 1-e^{-r^2/2}\n",
|
|
"$$\n",
|
|
"integrating over values $r > 0$ as it cannot be negative, and using a substitution with $z := t^2$.\n",
|
|
"\n",
|
|
"Now the inversion.\n",
|
|
"\n",
|
|
"$$\n",
|
|
"p := F_R(r) = 1-e^{-r^2/2} \\implies r = \\sqrt{-2\\ln{(1-p)}}\n",
|
|
"$$\n",
|
|
"for $p \\in [0, 1]$. Do note that we can also calculate $r = \\sqrt{-2\\ln{(p)}}$ as $p$ and $1 - p$ are identically distributed, also keeping the argument to $\\ln$ positive, saving just one calculation, although I do not use this in the following."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 12,
|
|
"id": "e4023f99",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "86173970c865da7b0cb8ab78ec4a87b6",
|
|
"grade": true,
|
|
"grade_id": "cell-9bf8873cce1d179c",
|
|
"locked": false,
|
|
"points": 15,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAgvElEQVR4nO3df3DcZbn38ffV0hqJKSK04NCWpp5qDSSsNQ0tOGURWgMiPYyP01IUxAOxKiAzwmPlIPMcZaSMjPJDzlN2tMKjrK0DFDNQAVHSyo9IU90TaLEa0yKh1JQApYWGEHo9f+yPbtL8+CbNZjff/bxmOt3vj3tz7VKuvXPvfV+3uTsiIhJe4/IdgIiI5JYSvYhIyCnRi4iEnBK9iEjIKdGLiITcEfkOoC/HHnusz5gxI99hiIiMGZs3b37V3Sf3da0gE/2MGTNoamrKdxgiImOGmb3Y3zUN3YiIhJwSvYhIyCnRi4iEXEGO0YtIOLz77ru0tbXR2dmZ71BCo6SkhKlTpzJhwoTAbZToRSRn2traKCsrY8aMGZhZvsMZ89ydjo4O2traKC8vD9xOQzcikjOdnZ0cc8wxSvIjxMw45phjhvwbkhK9iOSUkvzIGs77qUQvIhJygRK9mdWa2TYzazGzFQPcN9fM3jOz/zXUtiIiYZBIJFi/fv2Q20Wj0ZwtFB30y1gzGw/cCSwE2oBNZlbv7lv7uO9m4NGhthUZK2aseDjzeMfKz+YxEilUiUSCpqYmzj333HyHkhGkR18DtLh7q7t3AWuAxX3cdyVwP9A+jLYiIjnxy1/+kpqaGiKRCF/96lf505/+RFVVFZ2dnbz11lucdNJJPP/88zQ0NLBgwQIuuOACKioqWL58OQcOHADgscceY/78+cyZM4cvfOEL7Nu3D4BNmzZx2mmnccopp1BTU8OePXu44YYbWLt2LZFIhLVr1/LWW2/xla98hblz5/KJT3yC3/zmNwDs37+fpUuXUlVVxZIlS9i/f3/O3oMg0ytPAF7KOm4DTs2+wcxOAC4APg3MHUrbrOeoA+oApk+fHiAskdyLxWLE4/HM8a7WjoPXZr5MXV1dPsIak66++moSicSIPmckEuHWW2/t9/oLL7zA2rVreeqpp5gwYQJf//rX2bZtG+effz7XX389+/fv54tf/CInn3wyDQ0NPPvss2zdupUTTzyR2tpaHnjgAaLRKDfeeCOPP/44paWl3HzzzfzoRz9ixYoVLFmyhLVr1zJ37lzefPNNjjzySL73ve/R1NTET37yEwCuu+46Pv3pT7N69WreeOMNampqOPvss7nrrrs48sgjaW5uprm5mTlz5ozoe5MtSKLv6yve3hvN3gp8293f6/WNcJC2yZPuMSAGUF1drY1spSDE43ESiQSRSKTH+a727cTjcSX6Avf73/+ezZs3M3dusv+5f/9+pkyZwg033MDcuXMpKSnh9ttvz9xfU1PDzJkzAbjwwgt58sknKSkpYevWrZx++ukAdHV1MX/+fLZt28aHP/zhzHNPmjSpzxgee+wx6uvrueWWW4DklNN//vOfbNy4kauuugqAqqoqqqqqcvMmECzRtwHTso6nAjt73VMNrEkl+WOBc82sO2BbkYLRuwefTvINDQ3AwTH6XfEVJBIJotFo5t5ly5Yp8Q9goJ53rrg7l1xyCTfddFOP87t27WLfvn28++67dHZ2UlpaChw6ddHMcHcWLlzIr371qx7XmpubA011dHfuv/9+Pvaxjx1ybbSmngYZo98EzDKzcjObCCwF6rNvcPdyd5/h7jOA+4Cvu/uDQdqKFJJ0Dz4tEonwXEklM1Y83OOL2NKKaI9efiKR6PEBIYXhrLPO4r777qO9PfnV4WuvvcaLL75IXV0d3//+97nooov49re/nbn/2WefZfv27Rw4cIC1a9fyqU99innz5vHUU0/R0tICwNtvv83f/vY3Zs+ezc6dO9m0aRMAe/fupbu7m7KyMvbu3Zt5zs985jPccccduCcHKv7yl78AsGDBAu69914Ann/+eZqbm3P2Pgzao3f3bjO7guRsmvHAanffYmbLU9dXDbXtyIQukhvZPXjoOdMmrSxSyw5qM8edrZo5XIgqKiq48cYbWbRoEQcOHGDChAksXryYI444gmXLlvHee+9x2mmn8Yc//IFx48Yxf/58VqxYwXPPPZf5YnbcuHHcfffdXHjhhbzzzjsA3HjjjXz0ox9l7dq1XHnllezfv5/3v//9PP7445x55pmsXLmSSCTCd77zHb773e9y9dVXU1VVhbszY8YMHnroIb72ta9x6aWXUlVVRSQSoaamJmfvg6U/ZQpJdXW1a+MRyYdoNEpjawfHL1s5pHa74slEn26nqZdJL7zwAh//+MfzHUYgDQ0N3HLLLTz00EP5DmVQfb2vZrbZ3av7ul9FzaToZY/LJxIJmDRt4Ab96Grfnkn4mpEjhUQlEKToZY/LRyIRSiuiQ36O0oooE6ckqwmmZ+TI2BKNRsdEb3441KMXgT5n1gxFWaSWskhyzD7dqxcpFOrRi4iEnHr0IkBja8ewevIiY4F69FKUYrEY0WiUaDQ64svyATY2NlEyvZKS6ZXEYrERf36RoVCPXopSdmmD9KKokZL9Za5KJfQ00r81BZnGumPHDs477zyef/75HudvuOEGFixYwNlnn91nuwcffJCPfvSjVFRUjEis+aREL0XrcL+A7Y++mB0bvve97w14/cEHH+S8884bkUTf3d3NEUfkL91q6EaKVnpcXmPz4ffee+9x+eWXc9JJJ7Fo0SL279/Pl7/8Ze677z4AVqxYQUVFBVVVVVxzzTU8/fTT1NfXc+211xKJRPjHP/5BIpFg3rx5VFVVccEFF/D6668DyVLFVVVVzJ8/n2uvvZaTTz4ZgLvvvpsvfOELfO5zn2PRokXs27ePs846izlz5lBZWZkpV7xjxw5mz57NZZddxsknn8xFF13E448/zumnn86sWbN49tlnD/v1K9GLSOj9/e9/5xvf+AZbtmzhgx/8IPfff3/m2muvvca6devYsmULzc3NXH/99Zx22mmcf/75/PCHPySRSPCRj3yEiy++mJtvvpnm5mYqKyv5r//6LwAuvfRSVq1axTPPPMP48eN7/NxnnnmGe+65hz/84Q+UlJSwbt06/vznP/PEE0/wrW99K1P/pqWlhW9+85s0Nzfz17/+lXg8zpNPPsktt9zCD37wg8N+/Ur0IjmW/ZuDfnvIj/Ly8kwRuk9+8pPs2LEjc23SpEmUlJRw2WWX8cADD3DkkUce0n7Pnj288cYbnHHGGQBccsklbNy4kTfeeIO9e/dy2mmnAckKptkWLlzIhz70ISBZxfK6666jqqqKs88+m5dffpl//etfmfgqKysZN24cJ510EmeddRZmRmVlZY9Yh0uJXkRC733ve1/m8fjx4+nu7s4cH3HEETz77LN8/vOf58EHH6S2travp+jTYLXC0uWPAe699152797N5s2bSSQSHHfccXR2dh4S37hx4zLH48aN6xHrcCnRi0hR27dvH3v27OHcc8/l1ltvzUy3zS43fNRRR3H00Ufzxz/+EYBf/OIXnHHGGRx99NGUlZXR2NgIwJo1a/r9OXv27GHKlClMmDCBJ554ghdffDG3LyyLZt1IUYjFYlx14x2Z46727ZnaNLmWXewMirvgWSFW9dy7dy+LFy+ms7MTd+fHP/4xAEuXLuXyyy/n9ttv57777uOee+5h+fLlvP3228ycOZOf//znAPzsZz/j8ssvp7S0lGg0ylFHHdXnz7nooov43Oc+R3V1NZFIhNmzZ4/aa1SZYikK0WiUjY1NPZJ7aUU0Mw0yV/YmHuGtrQ2Z46727SyYV92j3n2YjaUyxcO1b98+PvCBDwCwcuVKXnnlFW677bac/kyVKRbpx8Qp5UOuM3+4sufUg+bVh9HDDz/MTTfdRHd3NyeeeCJ33313vkM6RKBEb2a1wG0kd4n6qbuv7HV9MfB94ADQDVzt7k+mru0A9gLvAd39feKIiIxFS5YsYcmSJfkOY0CDJnozGw/cCSwkudn3JjOrd/etWbf9Hqh3dzezKuDXQPYA1Jnu/uoIxi0iY4S7j9om2MVgOMPtQWbd1AAt7t7q7l3AGmBxrx+8zw/+9FKg8Ab+RWTUlZSU0NHRMazkJIdydzo6OigpKRlSuyBDNycAL2UdtwGn9r7JzC4AbgKmANlfrTvwmJk5cJe791nKz8zqgDqA6dOnBwpeRArb1KlTaWtrY/fu3fkOJTRKSkqYOnXqkNoESfR9/c51yMezu68D1pnZApLj9emScKe7+04zmwL8zsz+6u4b+2gfA2KQnHUT9AWISOGaMGEC5eWjM41V+hdk6KYNyN4teSqws7+bU0n8I2Z2bOp4Z+rvdmAdyaEgEREZJUES/SZglpmVm9lEYClQn32Dmf2bpb5tMbM5wESgw8xKzawsdb4UWAT0LAotkiO53lxkuBKJRCYubUoio2HQoRt37zazK4BHSU6vXO3uW8xseer6KuDzwMVm9i6wH1iSmoFzHMnhnPTPirv7Izl6LSI95HJzkeEqrYhS2fkcQObDp1hXycroCTSP3t3XA+t7nVuV9fhm4OY+2rUCpxxmjCLD1jlpGjvmXQtAWZ5jgeQCqh0kF1B1tmrxlIwOFTUTEQk5JXoRkZBTohcRCTklehGRkFOiFxEJOSV6EZGQU6IXEQk5JXoJlUJdDdsfrZKV0aBEL6GSXg0LEIlEKK2I5jWegZRWROmcNI3G1g42NjYRj8fzHZKElLYSlNAptNWw/cneZlBbDEouqUcvIhJySvQiIiGnRC8iEnJK9CIiIacvY2XMm7Hi4czjXa0deYxEpDCpRy8iEnJK9CIiIRco0ZtZrZltM7MWMztkwq+ZLTazZjNLmFmTmX0qaFsREcmtQcfozWw8cCewEGgDNplZvbtvzbrt90B9ap/YKuDXwOyAbUWGLRaLsSt+R+a4q307E6eU5zGi4dvY2ETJ9IP72t5+/ZXaT1ZGRJAefQ3Q4u6t7t4FrAEWZ9/g7vvc3VOHpYAHbStyOOLxOF3t2zPHE6eUF3TZg/6UVkR7fEB1tW9XSQQZMUFm3ZwAvJR13Aac2vsmM7sAuAmYAnx2KG1T7euAOoDp06cHCEskaeKUco5ftjLfYRyW7HIIoJIIMrKC9Oitj3N+yAn3de4+G/h34PtDaZtqH3P3anevnjx5coCwREQkiCCJvg2YlnU8FdjZ383uvhH4iJkdO9S2IiIy8oIk+k3ALDMrN7OJwFKgPvsGM/s3M7PU4znARKAjSFsREcmtQcfo3b3bzK4AHgXGA6vdfYuZLU9dXwV8HrjYzN4F9gNLUl/O9tk2R69FRET6EKgEgruvB9b3Orcq6/HNwM1B24qIyOjRylgRkZBTohcRCTlVr5QxKV2xUtUqRQanHr1IgWps7WDGiod7lGEWGQ716GXMya5vM5Zr2wymq317ZoVsbObLqnsjw6YevYw52fVtxmptm8Fk175R3Rs5XOrRy5gUhvo2A8mufaO6N3K41KMXEQk5JXoRkZBTohcRCTklehGRkFOiFxEJOSV6EZGQU6IXEQk5JXoRkZBTohcRCblAid7Mas1sm5m1mNkhy/TM7CIza079edrMTsm6tsPMnjOzhJk1jWTwIiIyuEFLIJjZeOBOYCHJzb43mVm9u2/Num07cIa7v25m5wAx4NSs62e6+6sjGLcUkVgs1qPWSyKRgEnT+m8QQolEgmg0CsCyZctU4EyGJEiPvgZocfdWd+8C1gCLs29w96fd/fXUYSMwdWTDlGIWj8eTyT0lEomEspBZf0orokQiESCZ8FXgTIYqSFGzE4CXso7b6Nlb7+0/gN9mHTvwmJk5cJe7x4YcpRS9SCRCQ0ND5riYarSXRWppWJksy5zu1YsMRZBEb32c8z5vNDuTZKL/VNbp0919p5lNAX5nZn919419tK0D6gCmT58eICwpJulNOERk6IIk+jYge0B0KrCz901mVgX8FDjH3TP7u7n7ztTf7Wa2juRQ0CGJPtXTjwFUV1f3+UEiUqyyt06cN/OYPEcjY02QMfpNwCwzKzezicBSoD77BjObDjwAfMnd/5Z1vtTMytKPgUXA8yMVvIiIDG7QHr27d5vZFcCjwHhgtbtvMbPlqeurgBuAY4D/NjOAbnevBo4D1qXOHQHE3f2RnLwSERHpU6Adptx9PbC+17lVWY8vAy7ro10rcErv8yIiMnq0MlZEJOSU6EVEQk6JXkQk5JToRURCToleClIsFiMajRKNRnuUPxDY2NhEyfTKzJ9YTIvNZWBK9FKQ4vE4GxubaGztoHPStKKqbTOQ0oooE6eUZ4672rer9o0MKtD0SpF8mDilnOOXrcx3GAWlLFJLWaQ2c7wrfkjVcJFDqEcvIhJySvQiIiGnRC8iEnJK9CIiIadELyISckr0IiIhp0QvIhJySvQiIiGnRC8yxqX309WeutIfJXopGLFYLFO/ZWNjU77DGTO62rezK76CXfEVqnsjfQqU6M2s1sy2mVmLmR2y5trMLjKz5tSfp83slKBtRdLi8Thd7duBZPkD1bcZXHbtG9W9kf4MWuvGzMYDdwILgTZgk5nVu/vWrNu2A2e4++tmdg4QA04N2FYkQ/Vthia79o3q3kh/gvToa4AWd2919y5gDbA4+wZ3f9rdX08dNgJTg7YVEZHcCpLoTwBeyjpuS53rz38Avx1qWzOrM7MmM2vavXt3gLBERCSIIIne+jjnfd5odibJRP/tobZ195i7V7t79eTJkwOEJSIiQQSpR98GTMs6ngrs7H2TmVUBPwXOcfeOobQVEZHcCdKj3wTMMrNyM5sILAXqs28ws+nAA8CX3P1vQ2krIiK5NWiP3t27zewK4FFgPLDa3beY2fLU9VXADcAxwH+bGUB3ahimz7Y5ei0iItKHQFsJuvt6YH2vc6uyHl8GXBa0rYiIjB6tjBURCTklehGRkAs0dCOSC7FYrMeS/UQiAZOm9d9ABpVIJIhGowAsW7aMurq6/AYkBUE9esmbeDyeTO4pkUhE9W0OQ2lFlEgkAiQTvureSJp69JJXkUiEhoaGzLFK7Q5fWaSWhpV3AGR69SKgRC95lq6lLiMj/V7uau1g3sxj8hyNFAoN3YiIhJwSvYhIyCnRi4iEnBK9iEjIKdGLiIScEr2ISMgp0YuIhJwSvYhIyGnBlIyq7Po2qm2TW9l1b0C1b4qZevQyqrLr26i2Te5k170B1b4pdurRy6jLrm+j8ge5URapZQe1mePO1hV5jEbyLVCP3sxqzWybmbWY2SH/Ysxstpk9Y2bvmNk1va7tMLPnzCxhZk0jFbiIiAQzaI/ezMYDdwILgTZgk5nVu/vWrNteA64C/r2fpznT3V89zFglJFTITGR0BenR1wAt7t7q7l3AGmBx9g3u3u7um4B3cxCjiIgchiCJ/gTgpazjttS5oBx4zMw2m1m/X/mbWZ2ZNZlZ0+7du4fw9CIiMpAgid76OOdD+Bmnu/sc4BzgG2a2oK+b3D3m7tXuXj158uQhPL2IiAwkSKJvA7InO08Fdgb9Ae6+M/V3O7CO5FCQiIiMkiCJfhMwy8zKzWwisBSoD/LkZlZqZmXpx8Ai4PnhBisiIkM36Kwbd+82syuAR4HxwGp332Jmy1PXV5nZ8UATMAk4YGZXAxXAscA6M0v/rLi7P5KTVyIFS6thC0P2Slmtki0ugRZMuft6YH2vc6uyHu8iOaTT25vAKYcToIx98XicjY1NTJxSDpOmaTVsHpRWRKnsfA4gszJZib54aGWsjIqJU8o5ftnKfIdRtMoitTSsvAOgR/0bKQ6qdSMiEnLq0YsUifRq5F2tHcybeUyeo5HRpB69iEjIKdGLiIScEr2ISMgp0YuIhJwSvYhIyGnWjYy47JWwoNWwhWhjYxMl0yszx7dff6UWUIWYevQy4rL3hQXtDVtoSiuiyVXKKV3t27WfbMipRy850TlpGjvmXZs5LstjLNJTWaSWssjB/WR3xbWfbNipRy8iEnJK9CIiIadELyISckr0IiIhp0QvIhJygRK9mdWa2TYzazGzQ76iN7PZZvaMmb1jZtcMpa2IiOTWoInezMYDdwLnkNwe8EIzq+h122vAVcAtw2grIRCLxYhGo0Sj0R5z6GVsSG8zGI1GicVi+Q5HRliQHn0N0OLure7eBawBFmff4O7t7r4JeHeobSUc0tsFNrZ20KntAseU0oookUgESCZ8LZ4KnyALpk4AXso6bgNODfj8h9NWxhhtFzg2lUVq2UFyAVVnq0ZXwyhIj976OOcBnz9wWzOrM7MmM2vavXt3wKcXEZHBBEn0bUB2RaqpwM6Azx+4rbvH3L3a3asnT54c8OlFRGQwQRL9JmCWmZWb2URgKVAf8PkPp62IiIyAQcfo3b3bzK4AHgXGA6vdfYuZLU9dX2VmxwNNwCTggJldDVS4+5t9tc3RaxERkT4Eql7p7uuB9b3Orcp6vIvksEygtiIiMnq0MlZEJOSU6GXYYrEYJdMrKZleycbGpnyHIyMkvftUyfRKLZ4KCSV6GbZ4PE5X+3YgOYdei6TGvuzdp7TzVHhohyk5LFokFS7Zu09p56nwUI9eRCTklOhFREJOiV5EJOSU6EVEQk5fxkpgsViMq268I3Pc1b49M0NDwildpz5t2bJl1NXV5S8gGRYlegksPZ0yndw1pTLcSiuivLW1gcbWDoDMVFol+rFHiV6GRNMpi0f2VEvQdMuxTGP0IiIhp0QvIhJySvQiIiGnRC8DisViRKNRotEoiUQi3+FInqng2dikRC8DisfjmQQfiUQ0y6aIqeDZ2KVZNzKozknT2DHvWgDK8hyL5I8Kno1dgXr0ZlZrZtvMrMXMDvkvbEm3p643m9mcrGs7zOw5M0uYmYqWi4iMskF79GY2HrgTWAi0AZvMrN7dt2bddg4wK/XnVOD/pv5OO9PdXx2xqEVEJLAgPfoaoMXdW929C1gDLO51z2Lg/3lSI/BBM/vwCMcqIiLDECTRnwC8lHXcljoX9B4HHjOzzWbW79ppM6szsyYza9q9e3eAsCRXtEWgBKEZOGNHkERvfZzzIdxzurvPITm88w0zW9DXD3H3mLtXu3v15MmTA4QluaItAmUwmoEztgSZddMGTMs6ngrsDHqPu6f/bjezdSSHgjYON2AZHappIwPRDJyxJUiPfhMwy8zKzWwisBSo73VPPXBxavbNPGCPu79iZqVmVgZgZqXAIuD5EYxfREQGMWiP3t27zewK4FFgPLDa3beY2fLU9VXAeuBcoAV4G7g01fw4YJ2ZpX9W3N0fGfFXIYclFov1+NU7kUjApGn9NxDpRXXrC5u59x5uz7/q6mpvatKXgKMlXd4gEolkzj1XUtmjRK1If/YmHqGy87nMcfrfUkNDQ/6CKkJmttndq/u6ppWxAvRc/QpaASvBlUVq2cHBTkFnq8bsC41q3YiIhJwSfZFSVUrJJc2xLyxK9EUqHo+zsbGJxtYOOidN01x5GTGaY194NEZfxDRXXnJBc+wLj3r0IiIhp0RfRDQuL/mQnmMfjUY1Xp8nGropIulx+YlTykHj8jIKSiuivLW1gcbWjkz9JC2kGn1K9EVG4/IymjReXxiU6ENMpQ2k0KhUQn4o0YdYemPvdGmDSCTCcyWV+Q1Kilb2MA4kp142tnbwg9bk1hU7Vn42n+GFmhJ9yMRiMa668Q4g+T/SxCnlKm0gBSF7GAeSQzld7dszQzqxmS+rd58jmnUTMto0RMYKLawaPerRh0D2WHwikdAXrjIm9P6iNnv8XmP3I0uJfgzq/SXrhg0bAHjftJM1bVLGpNKKaKbU8YYNG9iwYUOPf+NK/IdH9ejHiOzk3iOxp5RWRFU/XkJhb+IR3trakDl+56XkpnRnnHEGoKTfn4Hq0SvRF7DsL1bT/9jTyV2JXYpF9sYm6U6Okv6hDjvRm1ktcBvJrQR/6u4re1231PVzSW4l+GV3/3OQtn0ppkTfexgmW++eu5K7FLvs3n7vnn5vxfYhcFiJ3szGA38DFgJtJDcLv9Ddt2bdcy5wJclEfypwm7ufGqRtX8Zioh8oYQ+kr2GYbEruIn3rPcSTbbAPgYGM1Q+Iw03084H/4+6fSR1/B8Ddb8q65y6gwd1/lTreBkSBGYO17UtZWZl/8pOfDPLa8i69+KP30MpQKJmLjKyBPgQG0vv/43kzjxnJsHJqw4YNh7Vn7AnAS1nHbSR77YPdc0LAtgCYWR2Q/hjdt2HDhm0BYsulY4FXh9oo/Q9lqG1ee/QnQ243iob1XoSU3ouDQvtepP8/3vDSIDceVAjvxYn9XQiS6K2Pc71/DejvniBtkyfdY0DB1DA1s6b+Ph2Ljd6Lg/ReHKT34qBCfy+CJPo2ILsS1lRgZ8B7JgZoKyIiORSkBMImYJaZlZvZRGApUN/rnnrgYkuaB+xx91cCthURkRwatEfv7t1mdgXwKMkpkqvdfYuZLU9dXwWsJznjpoXk9MpLB2qbk1cy8gpmGKkA6L04SO/FQXovDiro96IgF0yJiMjIUfVKEZGQU6IXEQk5JfoAzOwaM3MzOzbfseSLmf3QzP5qZs1mts7MPpjvmEaTmdWa2TYzazGzot381MymmdkTZvaCmW0xs2/mO6Z8M7PxZvYXM3so37H0R4l+EGY2jWQJh3/mO5Y8+x1wsrtXkSxr8Z08xzNqUqU87gTOASqAC82sIr9R5U038C13/zgwD/hGEb8Xad8EXsh3EANRoh/cj4H/TT8LvYqFuz/m7t2pw0aSayKKRQ3Q4u6t7t4FrAEW5zmmvHD3V9IFC919L8kEd0J+o8ofM5sKfBb4ab5jGYgS/QDM7HzgZXf/n3zHUmC+Avw230GMov5KfBQ1M5sBfAL4U55DyadbSXYED+Q5jgEV/Q5TZvY4cHwfl/4TuA5YNLoR5c9A74W7/yZ1z3+S/PX93tGMLc8Cl/IoFmb2AeB+4Gp3fzPf8eSDmZ0HtLv7ZjOL5jmcARV9onf3s/s6b2aVQDnwP8ly+0wF/mxmNe6+axRDHDX9vRdpZnYJcB5wlhfXAowgZUCKhplNIJnk73X3B/IdTx6dDpyfKtNeAkwys1+6+xfzHNchtGAqIDPbAVS7e74r1OVFagOZHwFnuPvufMczmszsCJJfQJ8FvEyytMeyMbTKe8SkNhm6B3jN3a/OczgFI9Wjv8bdz8tzKH3SGL0E9ROgDPidmSXMbFW+AxotqS+h06U8XgB+XYxJPuV04EvAp1P/DhKpHq0UMPXoRURCTj16EZGQU6IXEQk5JXoRkZBTohcRCTklehGRkFOiFxEJOSV6EZGQ+/+bikqj+kgblwAAAABJRU5ErkJggg==\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAg7ElEQVR4nO3df3TcdZ3v8ee7v4zEBLW0i4e2NN2tWwMNY01DC552EKgBkS5nVUpwRVyIdUWXcy5cK+v23FWulHM5ihS8ddZl8a4M7S5QzIGuIGjo+iPSdJ0tpVCNabGhdlOC1BYaSuB9/5gffJNOkm+aTGbyndfjnJ7M9/v9fGY+M4T3fPL5fj7vj7k7IiISXZOK3QARESksBXoRkYhToBcRiTgFehGRiFOgFxGJuCnFbkA+p5xyis+dO7fYzRARmTC2b9/+orvPyHetJAP93LlzaW9vL3YzREQmDDN7frBrGroREYk4BXoRkYhToBcRibiSHKMXkWh4/fXX6erqore3t9hNiYyKigpmzZrF1KlTQ9dRoBeRgunq6qKqqoq5c+diZsVuzoTn7vT09NDV1UVNTU3oehq6EZGC6e3tZfr06QryY8TMmD59+oj/QlKgF5GCUpAfWyfyeSrQi4hEnAK9iMgYSqVSbNmyZcT14vF4wRaKhgr0ZtZoZrvNrMPM1gxRbrGZvWFmHxtpXZGJYO6aR3L/RPI50UBfSMMGejObDNwFXATUAleYWe0g5W4FHh1pXRGRQvn+979PQ0MDsViMz372s/zyl7+krq6O3t5eXnnlFc444wx27txJa2sry5Yt47LLLqO2tpbVq1fz5ptvAvDYY4+xdOlSFi1axMc//nGOHDkCwLZt2zjnnHM466yzaGho4NChQ6xdu5ZNmzYRi8XYtGkTr7zyCp/5zGdYvHgx73//+/nBD34AwNGjR1m1ahV1dXVcfvnlHD16tGCfQZjplQ1Ah7t3ApjZRmAlsGtAuS8ADwCLT6CuiETc9ddfTyqVGtPnjMVi3H777YNef/bZZ9m0aRM/+9nPmDp1Kn/zN3/D7t27ufTSS/nKV77C0aNH+eQnP8mZZ55Ja2srTz31FLt27eL000+nsbGRBx98kHg8zs0338zjjz9OZWUlt956K9/4xjdYs2YNl19+OZs2bWLx4sX88Y9/5KSTTuKrX/0q7e3t3HnnnQDcdNNNfOhDH+Luu+/m5ZdfpqGhgQsuuIDvfOc7nHTSSezYsYMdO3awaNGiMf1sgsIE+tOAfYHjLuDsYAEzOw24DPgQ/QP9sHVFSlkikSCZTOaOD3T2vHVt3gs0NzcXo1kS0hNPPMH27dtZvDgdlo4ePcrMmTNZu3YtixcvpqKigjvuuCNXvqGhgXnz5gFwxRVX8NOf/pSKigp27drFueeeC8CxY8dYunQpu3fv5j3veU/uuaurq/O24bHHHqOlpYXbbrsNSE85/d3vfsfWrVv54he/CEBdXR11dXWF+RAIF+jzzeUZuKP47cCX3P2NAVN/wtRNFzRrBpoB5syZE6JZIoWXTCZJpVLEYrF+54917yGZTCrQj8BQPe9CcXeuuuoqbrnlln7nDxw4wJEjR3j99dfp7e2lsrISOH7qopnh7lx44YXcd999/a7t2LEj1FRHd+eBBx7gz//8z4+7Nl5TT8PcjO0CZgeOZwH7B5SpBzaa2V7gY8C3zewvQtYFwN0T7l7v7vUzZuRNqSwyLoI3XNs6e4jFYrS2ttLa2sqpTes4tWkd02aGX5UoxXP++edz//33093dDcBLL73E888/T3NzM1/72te48sor+dKXvpQr/9RTT7Fnzx7efPNNNm3axAc/+EGWLFnCz372Mzo6OgB49dVX+fWvf82CBQvYv38/27ZtA+Dw4cP09fVRVVXF4cOHc8/54Q9/mPXr1+Oe7uP+6le/AmDZsmXce++9AOzcuZMdO3YU7HMI06PfBsw3sxrgBWAV0BQs4O6533ozuwd42N0fMrMpw9UVKSWJRIIDyfW542Pde2iDvLNsUqkU8Xg8d9zU1KQefompra3l5ptvZsWKFbz55ptMnTqVlStXMmXKFJqamnjjjTc455xz+PGPf8ykSZNYunQpa9as4emnn87dmJ00aRL33HMPV1xxBa+99hoAN998M+9973vZtGkTX/jCFzh69Chvf/vbefzxxznvvPNYt24dsViML3/5y/z93/89119/PXV1dbg7c+fO5eGHH+Zzn/scV199NXV1dcRiMRoaGgr2OVj2W2bIQmYXkx6emQzc7e7/28xWA7j7hgFl7yEd6O8frO5wr1dfX+/aeETGUzaQH0iu4Vj3nn499sraOFWxxn7lD6d+yCu7WnPHx7r3sGxJPa2trchbnn32Wd73vvcVuxmhtLa2ctttt/Hwww8XuynDyve5mtl2d6/PVz5UUjN33wJsGXBuwyBlPz1cXZFSNm1mDac2rRuyTFWssV/wP5DUEhEpXcpeKWUvOFwzsDc/EsGhHA3jTDzxeLzfUFyUKAWClL1kMsmx7j1AujdfWRsf8XNU1sZzM3NSqVS/KZkixaYevQjhhmuGUhVrpHVd+q+CqPYKZeJSj15EJOIU6EVEIk5DN1K2clMqA2kNxur5lsybPibPGTVjnfVz77qPDF9m714uueQSdu7c2e/82rVrWbZsGRdccEHeeg899BDvfe97qa2d+HkYFeilLI3VTJvBaAZO6fvqV7865PWHHnqISy65ZEwCfV9fH1OmFC/cauhGytJYzLQZTGVtnN7q2bR19rC1rV0zcErAG2+8wbXXXssZZ5zBihUrOHr0KJ/+9Ke5//77AVizZg21tbXU1dVxww038POf/5yWlhZuvPFGYrEYv/3tb0mlUixZsoS6ujouu+wy/vCHPwDpVMV1dXUsXbqUG2+8kTPPPBOAe+65h49//ON89KMfZcWKFRw5coTzzz+fRYsWsXDhwly64r1797JgwQKuueYazjzzTK688koef/xxzj33XObPn89TTz016vevHr2UrdHOtBlMcDGVFlKVht/85jfcd999/OM//iOf+MQneOCBB3LXXnrpJTZv3sxzzz2HmfHyyy/zzne+k0svvZRLLrmEj30svY9SXV0d69evZ/ny5axdu5Z/+Id/4Pbbb+fqq68mkUhwzjnnsGZN///ev/jFL9ixYwfvfve76evrY/PmzVRXV/Piiy+yZMkSLr30UgA6Ojr4t3/7NxKJBIsXLyaZTPLTn/6UlpYWvv71r/PQQw+N6v2rRy8ikVdTU5Nb5/CBD3yAvXv35q5VV1dTUVHBNddcw4MPPshJJ510XP1Dhw7x8ssvs3z5cgCuuuoqtm7dyssvv8zhw4c555xzgPQwXdCFF17Iu9/9biCdxfKmm26irq6OCy64gBdeeIH//u//zrVv4cKFTJo0iTPOOIPzzz8fM2PhwoX92nqi1KMXKbC2zp5+NyHD3ECUsfW2t70t93jy5Mn9dnOaMmUKTz31FE888QQbN27kzjvv5Mc//nGo5x0uV1g2/THAvffey8GDB9m+fTtTp05l7ty59Pb2Hte+SZMm5Y4nTZpEX19fqLYMRT16ESlrR44c4dChQ1x88cXcfvvtuV2wgumGTz75ZN71rnfxH//xHwD8y7/8C8uXL+dd73oXVVVVtLW1AbBx48ZBX+fQoUPMnDmTqVOn8pOf/ITnn3++sG8sQD16ERk3pfjXzOHDh1m5ciW9vb24O9/85jcBWLVqFddeey133HEH999/P9/73vdYvXo1r776KvPmzeOf//mfAfinf/onrr32WiorK4nH45x88sl5X+fKK6/kox/9KPX19cRiMRYsWDBu7zFUmuLxpjTFUmjxeJy2zp6C3IwNyt6MDb5OKQa7QplIaYpP1JEjR3jHO94BwLp16/j973/Pt771rYK+ZkHSFItMdAP3fk2lUlA9e/AKY+hY955+s2+012y0PPLII9xyyy309fVx+umnc8899xS7ScfRGL2Uhezer1mxWGxM584PprI23m8xVnavWYmOyy+/nFQqxc6dO3nkkUcoxa1Q1aOXspHd+zVrrJfj56MNStIzU8ZrE+xycCLD7aF69GbWaGa7zazDzI77TTWzlWa2w8xSZtZuZh8MXNtrZk9nr424hSIyYVVUVNDT03NCwUmO5+709PRQUVExonrD9ujNbDJwF3Ah0AVsM7MWd98VKPYE0OLubmZ1wL8CwVvK57n7iyNqmcgYGzifXQpv1qxZdHV1cfDgwWI3JTIqKiqYNWvWiOqEGbppADrcvRPAzDYCK4FcoHf3I4HylYC+vkWEqVOnUlMztgnjZOTCDN2cBuwLHHdlzvVjZpeZ2XPAI8BnApcceMzMtpvZoFMNzKw5M+zTrm9/EZGxEybQ57uLclyP3d03u/sC4C+ArwUunevui4CLgM+b2bJ8L+LuCXevd/f6UrxrLSIyUYUJ9F1AcMLxLGD/YIXdfSvwp2Z2SuZ4f+ZnN7CZ9FCQiIiMkzCBfhsw38xqzGwasApoCRYwsz+zzPwpM1sETAN6zKzSzKoy5yuBFUD/bV5ECiSRSBCPx4nH4/3m0BdbdlOSeDxOIpEodnOkDAx7M9bd+8zsOuBRYDJwt7s/Y2arM9c3AH8JfMrMXgeOApdnZuD8CbA58x0wBUi6+w8L9F5E+skukorFYsRiMZ6uWFjsJlFZG2dh79MAuS8frZKVQgu1YMrdtwBbBpzbEHh8K3BrnnqdwFmjbKPICQsukiqFqZVVsUb2kl5A1dtZfounpDi0MlYiTXPnRZTrRkQk8hToRUQiToFeRCTiFOhFiih7D0H3EaSQdDNWIicbNA909hS5JSKlQT16EZGIU49eIiWRSHAguR5I7+YU3N2pFAW3GdQWg1Io6tFLpCSTSY517wFg2syacdku8EQFtxnUFoNSSOrRS+RMm1nDqU3rit2MYQW3GSzHLQZl/KhHLyIScQr0IiIRp0AvIhJxCvQiIhGnQC8iEnEK9CIiERcq0JtZo5ntNrMOMztuHpiZrTSzHWaWMrN2M/tg2LoiIlJYw86jN7PJwF3AhaQ3Ct9mZi3uvitQ7AmgJbN9YB3wr8CCkHVFhOM3Sdm77iNFbI1ESZgFUw1AR2ZbQMxsI7ASyAVrdz8SKF8JeNi6IqORSCT6rShNpVJQPbt4DRqFYDoEUEoEGTthhm5OA/YFjrsy5/oxs8vM7DngEeAzI6mbqd+cGfZpP3jwYJi2i5BMJtna1k5bZw9tnT30Vs8u6bQHgwmmQwClRJCxFaZHb3nO+XEn3DcDm81sGfA14IKwdTP1E0ACoL6+Pm8ZkXwmSsqDoQTTIYBSIsjYCtOj7wKCfwvPAvYPVtjdtwJ/amanjLSuiIiMvTCBfhsw38xqzGwasApoCRYwsz8zM8s8XgRMA3rC1BURkcIadujG3fvM7DrgUWAycLe7P2NmqzPXNwB/CXzKzF4HjgKXu7sDeesW6L2IiEgeodIUu/sWYMuAcxsCj28Fbg1bV0RExo9WxoqIRJwCvYhIxCnQi4hEnAK9iEjEKdCLiEScAr1MOIlEgng8TjweT+e2iahUKpV7n4lEotjNkQlMgV4mnGQymQvwsVhsQua2GU5lbZxYLAakA77y3shohJpHL1JqYrEYra2tAP1S+0ZFVayR1nXrAYjH48VtjEx46tGLiEScevQyIQ3cpCOKsu/vQGcPS+ZNL3JrZCJTj15EJOIU6EVEIk6BXkQk4hToRUQiToFeRCTiFOhFRCJOgV5EJOJCzaM3s0bgW6S3A/yuu68bcP1K4EuZwyPA59z9vzLX9gKHgTeAPnevH5umS7lIJBL9UgCkUimonj14hQjK5r0BaGpqorm5ubgNkgll2B69mU0G7gIuAmqBK8ysdkCxPcByd68DvgYMzMB0nrvHFOTlRARz20B089sMRnlvZLTC9OgbgA537wQws43ASmBXtoC7/zxQvg2YNZaNFOmtns3eJTfmjquK2Jbxprw3MlphAv1pwL7AcRdw9hDl/xr498CxA4+ZmQPfcfe8+VbNrBloBpgzZ06IZomUD6VDkNEIE+gtzznPW9DsPNKB/oOB0+e6+34zmwn8yMyec/etxz1h+gsgAVBfX5/3+UVEZOTCzLrpAoJ3vmYB+wcWMrM64LvASnfvyZ539/2Zn93AZtJDQSIiMk7CBPptwHwzqzGzacAqoCVYwMzmAA8Cf+Xuvw6crzSzquxjYAWwc6waLyIiwxt26Mbd+8zsOuBR0tMr73b3Z8xsdeb6BmAtMB34tpnBW9Mo/wTYnDk3BUi6+w8L8k5ERCSvUPPo3X0LsGXAuQ2Bx9cA1+Sp1wmcNco2iojIKGhlrIhIxCnQi4hEnAK9iEjEac9YKUnB/DblmNtmKMG8N6DcNzI8BXopSclkkq1t7UybWQPVs8sqt81QKmvjLOx9OneczQGkQC9DUaCXkjVtZg2nNq0bvmAZqYo1spfG3HFv55oitkYmCo3Ri4hEnAK9iEjEKdCLiEScAr2ISMQp0IuIRJwCvYhIxCnQi4hEnAK9iEjEKdCLiEScAr2UjEQiQcWchVTMWcjWtvZiN2fCyOa+icfjJBKJYjdHSlCoQG9mjWa228w6zOy4NddmdqWZ7cj8+7mZnRW2rkhWMpnkWPceIJ3+QPlthldZGycWiwHpgJ9NBCcSNGyuGzObDNwFXEh6o/BtZtbi7rsCxfYAy939D2Z2EZAAzg5ZVyRH+W1GpirWSOu69QD9MlqKBIVJatYAdGS2BcTMNgIrgVywdvefB8q3AbPC1hWR0Zm75hEADnT2sGTe9CK3RkpRmKGb04B9geOuzLnB/DXw7yOta2bNZtZuZu0HDx4M0SwREQkjTKC3POc8b0Gz80gH+i+NtK67J9y93t3rZ8yYEaJZIiISRpihmy4guL3PLGD/wEJmVgd8F7jI3XtGUldERAonTI9+GzDfzGrMbBqwCmgJFjCzOcCDwF+5+69HUldERApr2B69u/eZ2XXAo8Bk4G53f8bMVmeubwDWAtOBb5sZQF9mGCZv3QK9FxERySPUVoLuvgXYMuDchsDja4BrwtYVEZHxo5WxIiIRp0AvIhJxoYZuRAohkUj0W7KfSqWgevbgFWRY2bw3AE1NTTQ3Nxe3QVIS1KOXokkmk+ngnhGLxZTfZhSU90YGox69FFUsFqO1tTV3nF3OLyOnvDcyGAV6Kaq2zh4F9zGkvDeSj4ZuREQiToFeRCTiFOhFRCJOgV5EJOIU6EVEIk6BXkQk4jS9UiSiBk5d3bvuI0VsjRSTAr2Mq2DaA6U8KKxj3Xs4kFyTO07Me0EpEcqUhm5kXAXTHijlQeFU1saZNrMmd3yse49SIpQx9ehl3AXTHmhVbGFUxRqpijXmjoM9eyk/oXr0ZtZoZrvNrMPMjvuNMbMFZvYLM3vNzG4YcG2vmT1tZikzax+rhouISDjD9ujNbDJwF3Ah6c2+t5lZi7vvChR7Cfgi8BeDPM157v7iKNsqIiInIEyPvgHocPdOdz8GbARWBgu4e7e7bwNeL0AbRURkFMKM0Z8G7AscdwFnj+A1HHjMzBz4jrsnRlBXIkgZK0XGV5hAb3nO+Qhe41x3329mM4Efmdlz7r71uBcxawaaAebMmTOCpxcRkaGEGbrpAoKTnWcB+8O+gLvvz/zsBjaTHgrKVy7h7vXuXj9jxoywTy8iIsMIE+i3AfPNrMbMpgGrgJYwT25mlWZWlX0MrAB2nmhjRURk5IYdunH3PjO7DngUmAzc7e7PmNnqzPUNZnYq0A5UA2+a2fVALXAKsNnMsq+VdPcfFuSdiIhIXqEWTLn7FmDLgHMbAo8PkB7SGeiPwFmjaaBMfEp7UBpSqVRuL9mmpialQygjSoEgBae0B8VXWRsnFosB6YCvdAjlRSkQZFz0Vs9m75IbAagqclvKUVWskdZ16wFyvXopH+rRi4hEnHr0ImUiu0jtQGcPS+ZNL3JrZDypRy8iEnEK9CIiEadALyIScQr0IiIRp5uxUhDB7JQHOnuK2BIRUaCXMZdIJDiQXJ87Pta9p9/+pVJ8wVWyoJWyUaehGxlzyWSSY917csfTZtZoNWwJCa6SBa2ULQfq0UtBTJtZw6lN64rdDMkjuEoWtFK2HCjQi5ShgfdQtIAq2jR0IyIScQr0IpLbx1d7+UaTAr2ISMQp0IuIRFyoQG9mjWa228w6zGxNnusLzOwXZvaamd0wkroiIlJYwwZ6M5sM3AVcRHof2CvMrHZAsZeALwK3nUBdiYBEIkHFnIVUzFnI1rb2YjdHRuhY9x4OJNdwILmGRCJR7ObIGAvTo28AOty9092PARuBlcEC7t7t7tuA10daV6IhuEhKC6QmlsraeG7l8rHuPVo8FUFh5tGfBuwLHHcBZ4d8/tB1zawZaAaYM2dOyKeXUqJFUhNTVayRqlgjAAeSGl2NojA9estzzkM+f+i67p5w93p3r58xY0bIpxcRkeGECfRdwOzA8Sxgf8jnH01dEREZA2EC/TZgvpnVmNk0YBXQEvL5R1NXRETGwLBj9O7eZ2bXAY8Ck4G73f0ZM1udub7BzE4F2oFq4E0zux6odfc/5qtboPciIiJ5hEpq5u5bgC0Dzm0IPD5AelgmVF0RERk/WhkrJyyRSBCPx4nH46RSqWI3R8ZIdlOSeDyuOfURoTTFcsKSySRb29rTc7CrZ2vufARU1sZ5ZVcrbZ09uXUR2nlq4lOgl1HR3Plo0Zz6aNLQjYhIxCnQi4hEnAK9iEjEKdCLiEScbsZKaIlEgi/evD53fKx7Ty7roURTdqplVlNTk2bhTEDq0UtowVTEoHTEUVdZGycWi+WOU6mUUhhPUOrRy4hoOmX5qIo1spfG3HFvp6ZbTlTq0YuIRJwCvYhIxCnQi4hEnAK9DEmJyyRICc8mJt2MlSEpcZlkKeHZxKVAL8PSTBsBJTybyEIN3ZhZo5ntNrMOMzvuv7Cl3ZG5vsPMFgWu7TWzp80sZWbtY9l4EREZ3rA9ejObDNwFXEh6s+9tZtbi7rsCxS4C5mf+nQ3838zPrPPc/cUxa7WIiIQWpkffAHS4e6e7HwM2AisHlFkJ/D9PawPeaWbvGeO2iojICQgT6E8D9gWOuzLnwpZx4DEz225mg965MbNmM2s3s/aDBw+GaJYUSiKRoGLOQirmLGRrm0bbJD/NwJk4wgR6y3POR1DmXHdfRHp45/Nmtizfi7h7wt3r3b1+xowZIZolhRLMaaN8NpJPMA+OcuCUvjCzbrqA2YHjWcD+sGXcPfuz28w2kx4K2nqiDZbxoZk2MpRgHhzlwCl9YXr024D5ZlZjZtOAVUDLgDItwKcys2+WAIfc/fdmVmlmVQBmVgmsAHaOYftFRGQYw/bo3b3PzK4DHgUmA3e7+zNmtjpzfQOwBbgY6ABeBa7OVP8TYLOZZV8r6e4/HPN3ISIigwq1YMrdt5AO5sFzGwKPHfh8nnqdwFmjbKMUWCKR6DfGmkqloHr24BVEBtja1k7FnIW54zu+8gWtmi0hynUjJJPJfnlsYrGYbsBKaJW18X47jR3r3qObsyVGKRAEgN7q2exdcmPuuKqIbZGJJZgaAZQeoRSpRy8iEnEK9GVKi6KkkLJj9hVzFmoxVQlQoC9TWhQlhRIcs9d4fWnQGH0Z06IoKQSlMy496tGXEe0WJcWgYZziU6AvI9ndoto6e+jVblEyDjSMUxo0dFNmNFwj40nDOKVBgT7CtOJVSk02tXFWU1OTVtCOAw3dRJhWvEopqayN01s9m7bOHto6e9ja1q6hnHGiHn3EBHvxqVRKK16lZORbQRvMkaP8OIWjHn3E6IarTBS6UTt+1KOPgIG9eN1wlYlg4I1a9e4LR4F+Ahp4k/XJJ58E4G2zzwT14mUCCv7OvrZvJ5/97Gf7/Y7rpu3oWDqVfGmpr6/39nblXwkKBvdsYF++fHnu+tMVC/uNf4pMVIdTP+SVXa2549f2pTely/6+K+jnZ2bb3b0+7zUF+tKVL7i/bfaZQLoHpMAu5SAY+LNBP/v/gYZ43jLqQG9mjcC3SG8l+F13XzfgumWuX0x6K8FPu/t/hqmbTzkF+oHDMEEK7iL95Qv6wb9sg8qt5z+qQG9mk4FfAxcCXaQ3C7/C3XcFylwMfIF0oD8b+Ja7nx2mbj4TMdAPFbAHauvsyT0e2EMZSMFdJL+BQzxBw30JDGWifkGMNtAvBf6Xu384c/xlAHe/JVDmO0Cru9+XOd4NxIG5w9XNp6qqyj/wgQ+EeW8FEwzGS+ZNH7bccAF7KArmImNrqC+BoYzmC6LYnnzyyUEDfZhZN6cB+wLHXaR77cOVOS1kXQDMrBnIfo0eefLJJ3eHaFshnQK8CPDkvmFKBmR/UUbitX07eenRO0dcbxzlPgvRZxEQ2c8iO2w6AqXwWZw+2IUwgd7ynBv4Z8BgZcLUTZ90TwAlk8PUzNoH+3YsN/os3qLP4i36LN5S6p9FmEDfBQQzYc0C9ocsMy1EXRERKaAwKRC2AfPNrMbMpgGrgJYBZVqAT1naEuCQu/8+ZF0RESmgYXv07t5nZtcBj5KeInm3uz9jZqsz1zcAW0jPuOkgPb3y6qHqFuSdjL2SGUYqAfos3qLP4i36LN5S0p9FSS6YEhGRsaPslSIiEadALyIScQr0IZjZDWbmZnZKsdtSLGb2f8zsOTPbYWabzeydxW7TeDKzRjPbbWYdZla2m5+a2Wwz+4mZPWtmz5jZ3xa7TcVmZpPN7Fdm9nCx2zIYBfphmNls0ikcflfsthTZj4Az3b2OdFqLLxe5PeMmk8rjLuAioBa4wsxqi9uqoukD/oe7vw9YAny+jD+LrL8Fni12I4aiQD+8bwL/k0EWepULd3/M3fsyh22k10SUiwagw9073f0YsBFYWeQ2FYW7/z6bsNDdD5MOcKcVt1XFY2azgI8A3y12W4aiQD8EM7sUeMHd/6vYbSkxnwH+vdiNGEeDpfgoa2Y2F3g/8MsiN6WYbifdEXyzyO0YUtnvMGVmjwOn5rn0d8BNwIrxbVHxDPVZuPsPMmX+jvSf7/eOZ9uKLHQqj3JhZu8AHgCud/c/Frs9xWBmlwDd7r7dzOJFbs6Qyj7Qu/sF+c6b2UKgBvivdLp9ZgH/aWYN7n5gHJs4bgb7LLLM7CrgEuB8L68FGGHSgJQNM5tKOsjf6+4PFrs9RXQucGkmTXsFUG1m33f3Txa5XcfRgqmQzGwvUO/uxc5QVxSZDWS+ASx394PFbs94MrMppG9Anw+8QDq1R9MEWuU9ZjKbDH0PeMndry9yc0pGpkd/g7tfUuSm5KUxegnrTqAK+JGZpcxsQ7EbNF4yN6GzqTyeBf61HIN8xrnAXwEfyvwepDI9Wilh6tGLiEScevQiIhGnQC8iEnEK9CIiEadALyIScQr0IiIRp0AvIhJxCvQiIhH3/wFF20gDtPdoFAAAAABJRU5ErkJggg==\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"def random_normal_pair():\n",
|
|
" '''Return two independent normal random variables.'''\n",
|
|
" phi = rng.random()*2*np.pi\n",
|
|
" r = inversion_sample(lambda p: np.sqrt(-2*np.log(1-p)))\n",
|
|
" x, y = r*np.cos(phi), r*np.sin(phi)\n",
|
|
" return x, y\n",
|
|
"\n",
|
|
"# Plotting\n",
|
|
"pdf = gaussian\n",
|
|
"samples = np.array([random_normal_pair() for _ in range(100000)])\n",
|
|
"for index in [0, 1]:\n",
|
|
" compare_plot(samples[:,index], pdf, -5, 5, 100) "
|
|
]
|
|
}
|
|
],
|
|
"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
|
|
}
|