979 lines
117 KiB
Plaintext
979 lines
117 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/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAXwklEQVR4nO3df3BV5Z3H8feXEBvB+KMIXYeoYAfFAPGWhixgF66ANLIq43YdILQqLaS0anWmdaWOZafKVBzdFmntItNS3NZIdhAtq2y1VJH6AyG0dyM/ik0BNaIlgiAgEQPf/SMRQsiPk+Te3OTJ5zWTMeec55775Yx87sNzn/Mcc3dERCQsPdJdgIiIJJ/CXUQkQAp3EZEAKdxFRAKkcBcRCVDPdL3xueee6wMGDEjX24uIdEkbN2583937ttQubeE+YMAAysrK0vX2IiJdkpm9GaWdhmVERAKkcBcRCZDCXUQkQGkbcxeRMHzyySdUVlZSXV2d7lKCkpWVRU5ODpmZmW16vcJdRNqlsrKS7OxsBgwYgJmlu5wguDt79uyhsrKSgQMHtukcGpYRkXaprq6mT58+CvYkMjP69OnTrn8NKdxFpN0U7MnX3muqcBcRCZDCXUSkDRKJBKtWrWr16+LxeIfcwKkvVNNs8eLFfGfeTyO3X3j3rRQXF6ewIhGJIpFIUFZWxqRJk9JdSqPUc0+zkpISjuzeEantkd07KCkpSXFFIl3Tb37zGwoKCojFYnzzm9/ktddeIy8vj+rqag4dOsSQIUPYtGkTa9asYcyYMVx33XXk5uYye/Zsjh07BsBzzz3HqFGjGD58ONdffz0HDx4EYMOGDYwePZrLLruMgoIC9u/fz9y5cyktLSUWi1FaWsqhQ4f4+te/zogRI/jCF77Ab3/7WwAOHz7M1KlTycvLY8qUKRw+fLhDrod67p3Aaf0G8g9F81ts917JnA6oRqTtbr/9dhKJRFLPGYvFWLBgQbNttm7dSmlpKS+//DKZmZl8+9vfZtu2bVx77bXcfffdHD58mK9+9asMHTqUNWvWsH79erZs2cKFF15IYWEhK1asIB6PM2/ePFavXk3v3r25//77+fGPf8ycOXOYMmUKpaWljBgxgg8//JBevXpxzz33UFZWxs9+9jMA7rrrLsaNG8eSJUvYt28fBQUFTJgwgUceeYRevXpRXl5OeXk5w4cPT+r1aYrCXUS6vD/84Q9s3LiRESNGALW95X79+jF37lxGjBhBVlYWCxcuPN6+oKCAiy66CIBp06bx0ksvkZWVxZYtW7j88ssBOHLkCKNGjWLbtm2cd955x8995plnNlrDc889x8qVK3nwwQeB2imib731FmvXruU73/kOAHl5eeTl5aXmIjSgcBeRpGmph50q7s6NN97Ifffdd9L+9957j4MHD/LJJ59QXV1N7969gVOnGZoZ7s6VV17J448/ftKx8vLySNMS3Z0nnniCSy655JRj6ZgqqjF3Eenyxo8fz/Lly9m9ezcAe/fu5c0336S4uJh7772X6dOnc+eddx5vv379enbs2MGxY8coLS3lS1/6EiNHjuTll1+moqICgI8++og33niDwYMHs2vXLjZs2ADAgQMHqKmpITs7mwMHDhw/55e//GV++tOf4u4A/PnPfwZgzJgxPPbYYwBs2rSJ8vLy1F8Q1HMXkQDk5uYyb948Jk6cyLFjx8jMzGTy5Mn07NmToqIijh49yujRo3n++efp0aMHo0aNYs6cObz++uvHv1zt0aMHS5cuZdq0aXz88ccAzJs3j4svvpjS0lJuvfVWDh8+zOmnn87q1au54oormD9/PrFYjO9///v84Ac/4PbbbycvLw93Z8CAATz99NN861vfYsaMGeTl5RGLxSgoKOiQa6JwF5EgTJkyhSlTpjR6LCMjg9deew2ANWvW0KtXL0pLS09pN27cuOM99PpGjBjBunXrTtnfsO0jjzxySpvTTz+dZcuWRfozJJOGZUREAqSeu4h0K/F4nHg8nu4yUk49dxGRACncRUQCpHAXEQmQwl1EJED6QlVEkmrAnGeSer6d8/+55TY7d3L11VezadOmk/bPnTuXMWPGMGHChEZf99RTT3HxxReTm5ublFo7E/XcRSRY99xzT5PBDrXhvmXLlqS8V01NTVLOkywKdxEJwtGjR5k1axZDhgxh4sSJHD58mJtuuonly5cDMGfOHHJzc8nLy+N73/ser7zyCitXruSOO+4gFovxt7/9jUQiwciRI8nLy+O6667jgw8+AGpvVsrLy2PUqFHccccdDB06FIClS5dy/fXXc8011zBx4kQOHjzI+PHjGT58OMOGDTu+7O/OnTsZPHgwM2fOZOjQoUyfPp3Vq1dz+eWXM2jQINavX5/069FiuJvZEjPbbWabmjg+3czK635eMbPLkl6liEgL/vrXv3LzzTezefNmzj77bJ544onjx/bu3cuTTz7J5s2bKS8v5+6772b06NFce+21PPDAAyQSCT7/+c9zww03cP/991NeXs6wYcP44Q9/CMCMGTNYtGgRr776KhkZGSe976uvvsqjjz7K888/T1ZWFk8++SR/+tOfeOGFF/jud797fK2ZiooKbrvtNsrLy/nLX/5CSUkJL730Eg8++CA/+tGPkn49ovTclwKFzRzfAYx19zzgXmBxEuoSEWmVgQMHEovFAPjiF7/Izp07jx8788wzycrKYubMmaxYsYJevXqd8vr9+/ezb98+xo4dC8CNN97I2rVr2bdvHwcOHGD06NEAFBUVnfS6K6+8ks9+9rNA7cqQd911F3l5eUyYMIF33nmHv//978frGzZsGD169GDIkCGMHz8eM2PYsGEn1ZosLYa7u68F9jZz/BV3/6Bucx2Qk6TaREQi+8xnPnP894yMjJPGwHv27Mn69ev5yle+wlNPPUVhYXP91ZN92vNuyqfLCAM89thjVFVVsXHjRhKJBJ/73Oeorq4+pb4ePXoc3+7Ro0dKxuuTPeb+DeB/mzpoZsVmVmZmZVVVVUl+axGRxh08eJD9+/czadIkFixYcPxpUfWX7T3rrLM455xz+OMf/wjAr3/9a8aOHcs555xDdnb28YXDmlsEbP/+/fTr14/MzExeeOEF3nzzzdT+wZqRtKmQZnYFteH+pabauPti6oZt8vPzm/84FJEuKcrUxY524MABJk+eTHV1Ne7OT37yEwCmTp3KrFmzWLhwIcuXL+fRRx9l9uzZfPTRR1x00UX86le/AuCXv/wls2bNonfv3sTjcc4666xG32f69Olcc8015OfnE4vFGDx4cIf9GRuylv7JAWBmA4Cn3X1oE8fzgCeBq9z9jShvnJ+f72VlZa0oNUzxeJx12/dEfobqyIv6sGbNmtQXJhLR1q1bufTSS9NdRkodPHiQM844A4D58+fz7rvv8tBDD6X8fRu7tma20d3zW3ptu3vuZnYBsAL4WtRgFxHpSp555hnuu+8+ampquPDCC1m6dGm6S2pRi+FuZo8DceBcM6sE/h3IBHD3RcBcoA/w87rnBNZE+VQREekqmnsQSGfVYri7+7QWjs8EZiatIhHpctw9LQ+BDlmUIfPm6A5VEWmXrKws9uzZ0+4wkhPcnT179pCVldXmc2jhMBFpl5ycHCorK9H05uTKysoiJ6fttw0p3EWkXTIzMxk4cGC6y5AGNCwjIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESA/r6GISiQTxeDxS26KiIoqLi1NbkIh0Sgr3LqR3bpxDW9awbvueFtse2b0DQOEu0k0p3LuQ7Fgh2bHCSG3fK5mT4mpEpDNrcczdzJaY2W4z29TEcTOzhWZWYWblZjY8+WWKiEhrRPlCdSnQXHfxKmBQ3U8x8J/tL0tERNqjxXB397XA3maaTAb+y2utA842s/OSVaCIiLReMqZC9gferrddWbfvFGZWbGZlZlZWVVWVhLcWEZHGJCPcrZF93lhDd1/s7vnunt+3b98kvLWIiDQmGeFeCZxfbzsH2JWE84qISBslI9xXAjfUzZoZCex393eTcF4REWmjFue5m9njQBw418wqgX8HMgHcfRGwCpgEVAAfATNSVayIiETTYri7+7QWjjtwc9IqEhGRdtPCYSIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gEKFK4m1mhmW0zswozm9PI8bPM7H/M7P/MbLOZzUh+qSIiElWL4W5mGcDDwFVALjDNzHIbNLsZ2OLulwFx4D/M7LQk1yoiIhH1jNCmAKhw9+0AZrYMmAxsqdfGgWwzM+AMYC9Qk+RapZUSiQTxeDxy+6KiIoqLi1NXkIh0mCjDMv2Bt+ttV9btq+9nwKXALuB14DZ3P9bwRGZWbGZlZlZWVVXVxpIlit65cWKxWOT2iUSCkpKS1BUkIh0qSs/dGtnnDba/DCSAccDngd+b2R/d/cOTXuS+GFgMkJ+f3/AckkTZsUJ2Uhi5ffX2U75KEZEuLErPvRI4v952DrU99PpmACu8VgWwAxicnBJFRKS1ooT7BmCQmQ2s+5J0KrCyQZu3gPEAZvY54BJgezILFRGR6FoclnH3GjO7BXgWyACWuPtmM5tdd3wRcC+w1Mxep3YY5053fz+FdYuISDOijLnj7quAVQ32Lar3+y5gYnJLExGRttIdqiIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAiPSBbuoe168rIumBYpLYL776V4uLiFFckIm2lnrsA0Ds3zmn9BkZqe2T3DkpKSlJckYi0R6Seu5kVAg8BGcAv3H1+I23iwAIgE3jf3ccmrUpJuexYIdmxwkht3yuZk+JqRKS9Wgx3M8sAHgauBCqBDWa20t231GtzNvBzoNDd3zKzfimqV0REIogyLFMAVLj7dnc/AiwDJjdoUwSscPe3ANx9d3LLFBGR1ogS7v2Bt+ttV9btq+9i4BwzW2NmG83shmQVKCIirRdlzN0a2eeNnOeLwHjgdOBVM1vn7m+cdCKzYqAY4IILLmh9tSIiEkmUnnslcH697RxgVyNtfufuh9z9fWAtcFnDE7n7YnfPd/f8vn37trVmERFpQZRw3wAMMrOBZnYaMBVY2aDNb4F/MrOeZtYL+Edga3JLFRGRqFoclnH3GjO7BXiW2qmQS9x9s5nNrju+yN23mtnvgHLgGLXTJTelsnAREWlapHnu7r4KWNVg36IG2w8ADySvNBERaSvdoSoiEiCFu4hIgBTuIiIBUriLiARI4S4iEiCFu4hIgBTuIiIBUriLiARI4S4iEiCFu4hIgBTuIiIBirS2jEhDiUSCeDweqW1RURHFxcWpLUhETqJwl1brnRvn0JY1rNu+p8W2R3bvAFC4i3Qwhbu0WnaskOxYYaS275XMSXE1ItIYjbmLiARI4S4iEiCFu4hIgBTuIiIBUriLiARI4S4iEiCFu4hIgBTuIiIBihTuZlZoZtvMrMLMmrwrxcxGmNlRM/vX5JUoIiKt1WK4m1kG8DBwFZALTDOz3Cba3Q88m+wiRUSkdaL03AuACnff7u5HgGXA5Eba3Qo8AexOYn0iItIGUcK9P/B2ve3Kun3HmVl/4DpgUfJKExGRtooS7tbIPm+wvQC4092PNnsis2IzKzOzsqqqqogliohIa0VZFbISOL/edg6wq0GbfGCZmQGcC0wysxp3f6p+I3dfDCwGyM/Pb/gBISIiSRIl3DcAg8xsIPAOMBUoqt/A3Qd++ruZLQWebhjs0n3pwR4iHa/FcHf3GjO7hdpZMBnAEnffbGaz645rnF2apAd7iKRHpId1uPsqYFWDfY2Gurvf1P6yJBR6sIdIeugOVRGRACncRUQCpHAXEQmQwl1EJEAKdxGRACncRUQCpHAXEQmQwl1EJEAKdxGRACncRUQCFGn5AZGOokXGRJJD4S6dhhYZE0kehbt0GlpkTCR5NOYuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAdI8d+mydDerSNMU7tIl6W5WkeYp3KVL0t2sIs2LNOZuZoVmts3MKszslL8pZjbdzMrrfl4xs8uSX6qIiETVYribWQbwMHAVkAtMM7PcBs12AGPdPQ+4F1ic7EJFRCS6KD33AqDC3be7+xFgGTC5fgN3f8XdP6jbXAfkJLdMERFpjSjh3h94u952Zd2+pnwD+N/GDphZsZmVmVlZVVVV9CpFRKRVooS7NbLPG21odgW14X5nY8fdfbG757t7ft++faNXKSIirRJltkwlcH697RxgV8NGZpYH/AK4yt1bnp8m0oFaMyceNC9eur4o4b4BGGRmA4F3gKlAUf0GZnYBsAL4mru/kfQqRdqhNXPiQfPiJQwthru715jZLcCzQAawxN03m9nsuuOLgLlAH+DnZgZQ4+75qStbJLrWzIkHzYuXMES6icndVwGrGuxbVO/3mcDM5JYmIiJtpYXDREQCpHAXEQmQ1pYRaYRWnJSuTuEu0oBWnJQQKNxFGtCKkxICjbmLiARI4S4iEiANy4i0k758lc5I4S7SDvryVTorhbtIO+jLV+msNOYuIhIg9dxFOpDG56WjKNxFOojG56UjKdxFOojG56UjKdxFOikN4Uh7KNxFOiEN4Uh7KdxFOqHWDuGoly8NKdxFujj18qUxCneRLi6VvXxQT7+rUriLdCOt6eUDfPz2Jl588UVKSkoitdcHQeehcBfpRlrTywc4kPhd5A8DfRB0Lgp3EWlSaz4M9EHQuUQKdzMrBB4CMoBfuPv8Bset7vgk4CPgJnf/U5JrFZFOrLN8EKRSV/qQaTHczSwDeBi4EqgENpjZSnffUq/ZVcCgup9/BP6z7r8iIqdI1QdBKnWmD5koovTcC4AKd98OYGbLgMlA/XCfDPyXuzuwzszONrPz3P3dpk66bdu2Vn1jH6pEIgFnnp/uMkQ6rdZ+T5AqneVDJqoo4d4feLvediWn9soba9MfOCnczawY+PTfNB+/+OKLm1pVbaj27z/3zfuvfj/dZXQS5wK6FrV0LU7QtTjhkiiNooS7NbLP29AGd18MLAYwszJ3z4/w/sHTtThB1+IEXYsTdC1OMLOyKO2iPKyjEqg/bpAD7GpDGxER6SBRwn0DMMjMBprZacBUYGWDNiuBG6zWSGB/c+PtIiKSWi0Oy7h7jZndAjxL7VTIJe6+2cxm1x1fBKyidhpkBbVTIWdEeO/Fba46PLoWJ+hanKBrcYKuxQmRroXVTnAREZGQ6AHZIiIBUriLiAQoreFuZteb2WYzO2Zm3XKak5kVmtk2M6sws2774EwzW2Jmu82s29/7YGbnm9kLZra17u/HbemuKV3MLMvM1pvZ/9Vdix+mu6Z0MrMMM/uzmT3dUtt099w3Af8CrE1zHWlRb2mHq4BcYJqZ5aa3qrRZCqT/NsTOoQb4rrtfCowEbu7G/198DIxz98uAGFBYNyOvu7oN2BqlYVrD3d23uvu2dNaQZseXdnD3I8CnSzt0O+6+Ftib7jo6A3d/99OF99z9ALV/mfunt6r08FoH6zYz63665SwQM8sB/hn4RZT26e65d3dNLdsgAoCZDQC+ALyW5lLSpm4oIgHsBn7v7t31WiwA/g04FqVxysPdzFab2aZGfrplD7WBSMs2SPdkZmcATwC3u/uH6a4nXdz9qLvHqL3zvcDMhqa5pA5nZlcDu919Y9TXpPxhHe4+IdXv0YVp2QZplJllUhvsj7n7inTX0xm4+z4zW0PtdzPd7Yv3y4FrzWwSkAWcaWa/cfevNvUCDcukV5SlHaSbqXv4zS+Bre7+43TXk05m1tfMzq77/XRgAvCXtBaVBu7+fXfPcfcB1ObE880FO6R/KuR1ZlYJjAKeMbNn01lPR3P3GuDTpR22Av/t7pvTW1V6mNnjwKvAJWZWaWbfSHdNaXQ58DVgnJkl6n4mpbuoNDkPeMHMyqntDP3e3VucBihafkBEJEgalhERCZDCXUQkQAp3EZEAKdxFRAKkcBcRCZDCXUQkQAp3EZEA/T/OJ56bDBK4PQAAAABJRU5ErkJggg==\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/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAbbUlEQVR4nO3de3BV5b3/8fcXCI0gXspFHWMldrAaIKQ0pFwciAo0elDGaR1uVmsPplgv9fxaf0XH4tRyKo5Oq6jnIKMUT2swZ1QoY1GpVcQbQrDbCCiWAtYUbVIQDEjEwPf8sTd0G5LslWRt9k7W5zWzh73X86yV734G+GTdnmXujoiIRFe3TBcgIiKZpSAQEYk4BYGISMQpCEREIk5BICIScT0yXUBz+vXr5wMHDsx0GSIincb69ev/6e7927NuVgbBwIEDqaqqynQZIiKdhpm93951dWhIRCTiFAQiIhGnIBARibisPEcgIl3D559/Tk1NDQ0NDZkupcvIzc0lLy+PnJyc0LapIBCRtKmpqaFPnz4MHDgQM8t0OZ2eu7Nz505qamrIz88Pbbs6NCQiadPQ0EDfvn0VAiExM/r27Rv6HpaCQETSSiEQrnSMp4JARCTiFAQiIiGKxWKsWLGizeuVlpZm7EZanSyOoIULF3Lj3PsD959/2w2Ul5ensSKRriMWi1FVVcXFF1+c6VIC0x5BBFVUVHCgdlugvgdqt1FRUZHmikTS53e/+x0lJSUUFRXxgx/8gDfeeIPCwkIaGhrYt28fgwcPZsOGDaxatYqxY8dy2WWXUVBQwKxZszh06BAAK1euZNSoUQwfPpzLL7+cvXv3ArBu3TpGjx7NsGHDKCkpYc+ePcyZM4fKykqKioqorKxk3759fP/732fEiBF8/etf5/e//z0A+/fvZ+rUqRQWFjJlyhT279+fsTHSHkFE9RyQz6nT56Xs91HF7GNQjUTBTTfdRCwWC3WbRUVF3HvvvS22v/POO1RWVvLqq6+Sk5PDD3/4QzZv3syll17Kbbfdxv79+7niiisYMmQIq1atYu3atWzatIkzzzyTsrIynnrqKUpLS5k7dy7PP/88vXv35q677uJXv/oVs2fPZsqUKVRWVjJixAg++eQTevXqxR133EFVVRUPPPAAALfeeisXXHABixYtYvfu3ZSUlDB+/HgeeughevXqRXV1NdXV1QwfPjzUsWkLBYGIdFl/+tOfWL9+PSNGjADiv4UPGDCAOXPmMGLECHJzc5k/f/6R/iUlJZx11lkATJs2jVdeeYXc3Fw2bdrEmDFjADhw4ACjRo1i8+bNnHbaaUe2fcIJJzRbw8qVK1m+fDn33HMPEL+k9m9/+xurV6/mxhtvBKCwsJDCwsL0DEIACgIROSZa+809Xdydq666ijvvvPMLyz/66CP27t3L559/TkNDA7179waOvjTTzHB3JkyYwJIlS77QVl1dHehSTnfnySef5Gtf+9pRbdlyaa3OEYhIl3XhhRfyxBNPUFtbC8CuXbt4//33KS8v5xe/+AUzZszgpz/96ZH+a9euZdu2bRw6dIjKykrOO+88Ro4cyauvvsqWLVsA+PTTT3nvvfc455xz2LFjB+vWrQOgvr6exsZG+vTpQ319/ZFtfutb3+L+++/H3QH485//DMDYsWN57LHHANiwYQPV1dXpH5AWaI9ARLqsgoIC5s6dy8SJEzl06BA5OTlMnjyZHj16MH36dA4ePMjo0aN54YUX6NatG6NGjWL27Nm8/fbbR04cd+vWjcWLFzNt2jQ+++wzAObOncvZZ59NZWUlN9xwA/v37+e4447j+eef5/zzz2fevHkUFRVxyy238LOf/YybbrqJwsJC3J2BAwfy9NNPc+2113L11VdTWFhIUVERJSUlGRsnBYGIdGlTpkxhypQpzbZ1796dN954A4BVq1bRq1cvKisrj+p3wQUXHPnNP9mIESNYs2bNUcub9n3ooYeO6nPcccfx+OOPB/oO6ZYyCMxsETAJqHX3Ic203wzMSNreuUB/d99lZtuBeuAg0OjuxWEVLiIi4QhyjmAxUNZSo7vf7e5F7l4E3AK85O67krqcn2hXCIhI1iotLeXpp5/OdBkZkTII3H01sCtVv4RpwJKUvUREJGuEdtWQmfUivufwZNJiB1aa2Xoza3WOAjMrN7MqM6uqq6sLqywREUkhzMtHLwFebXJYaIy7DwcuAq4zs7EtrezuC9292N2L+/fvH2JZIiLSmjCDYCpNDgu5+47En7XAUiBz10eJiEizQrl81MxOBMYBVyQt6w10c/f6xPuJwB1h/DwR6ZwGzv5DqNvbPu/fUvfZvp1JkyaxYcOGLyyfM2cOY8eOZfz48c2ut2zZMs4++2wKCgpCqTWbBbl8dAlQCvQzsxrgdiAHwN0XJLpdBqx0931Jq54CLE3cQt0DqHD3Z8MrXUSk/e64o/XfS5ctW8akSZNCCYLGxkZ69Mje27aCXDU0zd1Pc/ccd89z90fcfUFSCODui919apP1trr7sMRrsLv/Zzq+gIhIKgcPHuSaa65h8ODBTJw4kf379/O9732PJ554AoDZs2dTUFBAYWEhP/nJT3jttddYvnw5N998M0VFRfz1r38lFosxcuRICgsLueyyy/j444+B+M1jhYWFjBo1iptvvpkhQ+K3Wy1evJjLL7+cSy65hIkTJ7J3714uvPBChg8fztChQ49MR719+3bOOeccZs6cyZAhQ5gxYwbPP/88Y8aMYdCgQaxduzbt46O5hkSky/vLX/7Cddddx8aNGznppJN48sl/Xdy4a9culi5dysaNG6murua2225j9OjRXHrppdx9993EYjG++tWvcuWVV3LXXXdRXV3N0KFD+fnPfw7A1VdfzYIFC3j99dfp3r37F37u66+/zqOPPsoLL7xAbm4uS5cu5c033+TFF1/kxz/+8ZH5h7Zs2cKPfvQjqqureffdd6moqOCVV17hnnvu4Ze//GXax0dBICJdXn5+PkVFRQB84xvfYPv27UfaTjjhBHJzc5k5cyZPPfUUvXr1Omr9PXv2sHv3bsaNGwfAVVddxerVq9m9ezf19fWMHj0agOnTp39hvQkTJvDlL38ZiM9Ceuutt1JYWMj48eP5+9//zj/+8Y8j9Q0dOpRu3boxePBgLrzwQsyMoUOHfqHWdFEQiEiX96UvfenI++7du9PY2Hjkc48ePVi7di3f/va3WbZsGWVlLU6kcJTDv9G35PD01gCPPfYYdXV1rF+/nlgsximnnEJDQ8NR9XXr1u3I527dun2h1nRREIhIpO3du5c9e/Zw8cUXc++99x55ilrydNInnngiJ598Mi+//DIAv/3tbxk3bhwnn3wyffr0OTLxXGuTyO3Zs4cBAwaQk5PDiy++yPvvv5/eL9YG2XsaW0S6nCCXex5r9fX1TJ48mYaGBtydX//61wBMnTqVa665hvnz5/PEE0/w6KOPMmvWLD799FPOOussfvOb3wDwyCOPcM0119C7d29KS0s58cQTm/05M2bM4JJLLqG4uJiioiLOOeecY/YdU7FUuzaZUFxc7FVVVZkuo8sqLS1lzdadgZ9ZPPKsvqxatSr9hUmX884773Duuedmuoy02rt3L8cffzwA8+bN48MPP+S+++5L689sblzNbH17J/fUHoGISAf84Q9/4M4776SxsZEzzzyTxYsXZ7qkNlMQiIh0QGsPvuksdLJYRNIqGw8/d2bpGE8FgYikTW5uLjt37lQYhMTd2blzJ7m5uaFuV4eGRCRt8vLyqKmpQc8YCU9ubi55eXmhblNBICJpk5OTQ35+fqbLkBR0aEhEJOIUBCIiEacgEBGJOAWBiEjEKQhERCJOQSAiEnEKAhGRiEsZBGa2yMxqzWxDC+2lZrbHzGKJ15yktjIz22xmW8xsdpiFi4hIOILsESwGUj2y52V3L0q87gAws+7Ag8BFQAEwzcwKOlKsiIiEL2UQuPtqYFc7tl0CbHH3re5+AHgcmNyO7YiISBqFdY5glJm9ZWbPmNngxLLTgQ+S+tQkljXLzMrNrMrMqjQviYjIsRNGELwJnOnuw4D7gWWJ5dZM3xanIHT3he5e7O7F/fv3D6EsEREJosNB4O6fuPvexPsVQI6Z9SO+B3BGUtc8YEdHf56IiISrw0FgZqeamSXelyS2uRNYBwwys3wz6wlMBZZ39OeJiEi4Uk5DbWZLgFKgn5nVALcDOQDuvgD4DnCtmTUC+4GpHn8KRaOZXQ88B3QHFrn7xrR8CxERabeUQeDu01K0PwA80ELbCmBF+0oTEZFjQXcWi4hEnIJARCTiFAQiIhGnIBARiTgFgYhIxCkIREQiTkEgIhJxCgIRkYhTEIiIRJyCQEQk4hQEIiIRpyAQEYk4BYGISMSlnH1UOoeFCxdy49z7A/U9ULuNngPy01yRiHQW2iPoIioqKjhQuy1Q354D8uldUJregkSk09AeQRfSc0A+p06fl+kyRKST0R6BiEjEKQhERCJOQSAiEnEpg8DMFplZrZltaKF9hplVJ16vmdmwpLbtZva2mcXMrCrMwkVEJBxB9ggWA2WttG8Dxrl7IfALYGGT9vPdvcjdi9tXooiIpFPKq4bcfbWZDWyl/bWkj2uAvBDqEhGRYyTscwT/DjyT9NmBlWa23szKW1vRzMrNrMrMqurq6kIuS0REWhLafQRmdj7xIDgvafEYd99hZgOAP5rZu+6+urn13X0hicNKxcXFHlZdIiLSulD2CMysEHgYmOzuOw8vd/cdiT9rgaVASRg/T0REwtPhIDCzrwBPAd919/eSlvc2sz6H3wMTgWavPBIRkcxJeWjIzJYApUA/M6sBbgdyANx9ATAH6Av8l5kBNCauEDoFWJpY1gOocPdn0/AdRESkA4JcNTQtRftMYGYzy7cCw45eQ0REsonuLBYRiTjNPiopxWIxSktLA/WdPn065eWtXiksIllGQSCt6l1QytCGtwP1jcViAAoCkU5GQSCt6lNUxvZWZxj5l4ats9NcjYikg84RiIhEnIJARCTiFAQiIhGnIBARiTgFgYhIxCkIREQiTkEgIhJxCgIRkYhTEIiIRJyCQEQk4hQEIiIRpyAQEYk4BYGISMQpCEREIk5BICIScSmDwMwWmVmtmW1ood3MbL6ZbTGzajMbntRWZmabE22arF5EJAsF2SNYDK0+meQiYFDiVQ78N4CZdQceTLQXANPMrKAjxYqISPhSBoG7rwZ2tdJlMvA/HrcGOMnMTgNKgC3uvtXdDwCPJ/qKiEgWCeMcwenAB0mfaxLLWlreLDMrN7MqM6uqq6sLoSwREQkijCCwZpZ5K8ub5e4L3b3Y3Yv79+8fQlkiIhJEGA+vrwHOSPqcB+wAerawXEREskgYewTLgSsTVw+NBPa4+4fAOmCQmeWbWU9gaqKviIhkkZR7BGa2BCgF+plZDXA7kAPg7guAFcDFwBbgU+DqRFujmV0PPAd0Bxa5+8Y0fAcREemAlEHg7tNStDtwXQttK4gHhYiIZCndWSwiEnEKAhGRiFMQiIhEnIJARCTiFAQiIhGnIBARiTgFgYhIxCkIREQiTkEgIhJxCgIRkYhTEIiIRJyCQEQk4hQEIiIRpyAQEYk4BYGISMQpCEREIk5BICIScWE8vF7kiFgsRmlpaaC+06dPp7y8PL0FiUhKCgIJTe+CUvZtWsWarTtT9j1Quw1AQSCSBQIFgZmVAfcRfwj9w+4+r0n7zcCMpG2eC/R3911mth2oBw4Cje5eHFLtkmX6FJXRp6gsUN+PKmanuRoRCSplEJhZd+BBYAJQA6wzs+XuvulwH3e/G7g70f8S4D/cfVfSZs5393+GWrmIiIQiyMniEmCLu2919wPA48DkVvpPA5aEUZyIiKRfkCA4Hfgg6XNNYtlRzKwXUAY8mbTYgZVmtt7MWjwgbGblZlZlZlV1dXUByhIRkTAECQJrZpm30PcS4NUmh4XGuPtw4CLgOjMb29yK7r7Q3Yvdvbh///4ByhIRkTAECYIa4Iykz3nAjhb6TqXJYSF335H4sxZYSvxQk4iIZIkgQbAOGGRm+WbWk/h/9subdjKzE4FxwO+TlvU2sz6H3wMTgQ1hFC4iIuFIedWQuzea2fXAc8QvH13k7hvNbFaifUGi62XASnffl7T6KcBSMzv8syrc/dkwv4CIiHRMoPsI3H0FsKLJsgVNPi8GFjdZthUY1qEKRUQkrTTXkIhIxCkIREQiTkEgIhJxCgIRkYhTEIiIRJyCQEQk4hQEIiIRpyAQEYk4BYGISMQpCEREIk5BICIScQoCEZGIUxCIiEScgkBEJOIUBCIiEacgEBGJuEAPphFJh9Vrqsj9ytBAfeffdgPl5eVprkgkmrRHIBnRu6CUngPyA/U9ULuNioqKNFckEl3aI5CM6FNURp+iskB9P6qYneZqRKIt0B6BmZWZ2WYz22JmR/2rNLNSM9tjZrHEa07QdUVEJLNS7hGYWXfgQWACUAOsM7Pl7r6pSdeX3X1SO9cVEZEMCbJHUAJscfet7n4AeByYHHD7HVlXRESOgSBBcDrwQdLnmsSypkaZ2Vtm9oyZDW7juphZuZlVmVlVXV1dgLJERCQMQYLAmlnmTT6/CZzp7sOA+4FlbVg3vtB9obsXu3tx//79A5QlIiJhCBIENcAZSZ/zgB3JHdz9E3ffm3i/Asgxs35B1hURkcwKEgTrgEFmlm9mPYGpwPLkDmZ2qplZ4n1JYrs7g6wrIiKZlfKqIXdvNLPrgeeA7sAid99oZrMS7QuA7wDXmlkjsB+Y6u4ONLtumr6LiIi0Q6AbyhKHe1Y0WbYg6f0DwANB1xURkeyhKSZERCJOQSAiEnEKAhGRiFMQiIhEnIJARCTiNA21dAqxWIzS0tLA/adPn64H2YgEpCCQrNe7oJR9m1axZuvOQP0P1G4DUBCIBKQgkKzXlofYgB5kI9JWOkcgIhJxCgIRkYhTEIiIRJyCQEQk4hQEIiIRpyAQEYk4BYGISMQpCEREIk5BICIScbqzWLqktsxNpHmJJOoUBNLltGVuIs1LJBIwCMysDLiP+APoH3b3eU3aZwA/TXzcC1zr7m8l2rYD9cBBoNHdi8MpXaR5bZmbSPMSiQQIAjPrDjwITABqgHVmttzdNyV12waMc/ePzewiYCHwzaT28939nyHWLSIiIQlysrgE2OLuW939APA4MDm5g7u/5u4fJz6uAfLCLVNERNIlSBCcDnyQ9Lkmsawl/w48k/TZgZVmtt7MWjwQa2blZlZlZlV1dXUByhIRkTAEOUdgzSzzZjuanU88CM5LWjzG3XeY2QDgj2b2rruvPmqD7guJH1KiuLi42e2LiEj4guwR1ABnJH3OA3Y07WRmhcDDwGR3P3K5hrvvSPxZCywlfqhJRESyRJA9gnXAIDPLB/4OTAWmJ3cws68ATwHfdff3kpb3Brq5e33i/UTgjrCKFwmD7jmQqEsZBO7eaGbXA88Rv3x0kbtvNLNZifYFwBygL/BfZgb/ukz0FGBpYlkPoMLdn03LNxFpB91zIBLwPgJ3XwGsaLJsQdL7mcDMZtbbCgzrYI0iaaN7DkQ015CISOQpCEREIk5zDYm0gU4sS1ekIBAJSCeWpatSEIgEpBPL0lXpHIGISMRpj0AkTXQ+QToLBYFIGuh8gnQmCgKRNND5BOlMFAQiWUCHkSSTFAQiGda7oJShDW8H6huLxQAdRpJwKQhEMqxPURnbCXYYqWGrDiNJ+BQEIp1MWw4jgQ4lSWoKApFOpC1XIwF89sEGXnrpJSoqKgL1V2hEk4JApBNpy9VIAPWxZ3X+QVJSEIh0YW09/6Crl6JJQSAiQNsOO+mQU9eiIBARoG2Hnepjz7Y5NG6ce3+gbY88q2+gfqCACYuCQETarD2hEVS6ToS3VZRCJlAQmFkZcB/xh9c/7O7zmrRbov1i4FPge+7+ZpB1RaRra+sJ7qDaslfSVukOmWyTMgjMrDvwIDABqAHWmdlyd9+U1O0iYFDi9U3gv4FvBlxXRKTN0hUwkN6QyUZB9ghKgC3uvhXAzB4HJgPJ/5lPBv7H3R1YY2YnmdlpwMAA6x5l8+bNbbphRhKX/p1wRqbLEOkS0hky6fL+XZPavW6QIDgd+CDpcw3x3/pT9Tk94LoAmFk5cPiA3GcvvfTShgC1ZVI/4J+ZLuIL9uxp7i9D9tXZPNUZLtUZrs5Q59fau2KQILBmlnnAPkHWjS90XwgsBDCzKncvDlBbxnSGGkF1hk11hkt1hsfMqtq7bpAgqAGSjznkATsC9ukZYF0REcmgIM8sXgcMMrN8M+sJTAWWN+mzHLjS4kYCe9z9w4DriohIBqXcI3D3RjO7HniO+CWgi9x9o5nNSrQvAFYQv3R0C/HLR69ubd0AdS1sz5c5xjpDjaA6w6Y6w6U6w9PuGi1+oY+IiERVkENDIiLShSkIREQiLmNBYGZlZrbZzLaY2VHP30uceJ6faK82s+FZWmepme0xs1jiNScDNS4ys1oza/beiyway1R1ZnwsE3WcYWYvmtk7ZrbRzH7UTJ+Mj2nAOjM6pmaWa2ZrzeytRI0/b6ZPNoxlkDqz4u9nopbuZvZnM3u6mba2j6e7H/MX8RPHfwXOIn6J6VtAQZM+FwPPEL8XYSTwRpbWWQo8nYlxTKphLDAc2NBCe8bHMmCdGR/LRB2nAcMT7/sA72Xp388gdWZ0TBPjc3zifQ7wBjAyC8cySJ1Z8fczUcv/Ayqaq6c945mpPYIj01a4+wHg8NQTyY5MW+Hua4DD01ZkW50Z5+6rgV2tdMmGsQxSZ1Zw9w89MWmiu9cD7xC/Sz5Zxsc0YJ0ZlRifvYmPOYlX0ytUsmEsg9SZFcwsD/g34OEWurR5PDMVBC1NSdHWPukWtIZRiV3KZ8xs8LEprU2yYSyDyqqxNLOBwNeJ/4aYLKvGtJU6IcNjmjiMEQNqgT+6e1aOZYA6ITv+ft4L/H/gUAvtbR7PTAVBR6atOJaC1PAmcKa7DwPuB5alu6h2yIaxDCKrxtLMjgeeBG5y90+aNjezSkbGNEWdGR9Tdz/o7kXEZxYoMbMhTbpkxVgGqDPjY2lmk4Bad1/fWrdmlrU6npkKgo5MW3EspazB3T85vEvp7iuAHDPrd+xKDCQbxjKlbBpLM8sh/p/rY+7+VDNdsmJMU9WZTWPq7ruBVXDUQ5SzYiwPa6nOLBnLMcClZrad+KHqC8zsd036tHk8MxUEHZm2IqvqNLNTzcwS70uIj2m2TWKeDWOZUraMZaKGR4B33P1XLXTL+JgGqTPTY2pm/c3spMT744DxwLtNumXDWKasM9NjCeDut7h7nrsPJP7/0QvufkWTbm0ez4w8qtI7MG1FFtb5HeBaM2sE9gNTPXHq/lgxsyXEr2joZ2Y1wO3ET3ZlzVgGrDPjY5kwBvgu8HbimDHArcBXkmrNhjENUmemx/Q04FGLP6SqG/C/7v50tv1bD1hnpseyRR0dT00xISIScbqzWEQk4hQEIiIRpyAQEYk4BYGISMQpCEREIk5BICIScQoCEZGI+z+aMCvPDBrTTgAAAABJRU5ErkJggg==\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/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAa1ElEQVR4nO3df3BV5b3v8feHCKYiqMViHYIFPViMEHNpjCgOxl8UrZrj9HT41WrpkQxttTq39RYdS1v1tjJlWktrD820qKeVwhkVy1haf1Qd/MUhoc0JPxSbItYUWxArQg1i8Hv/yIa7CTtkJSTsDevzmsmw1nqeZ+3vfmDyYT1777UVEZiZWXr1yXcBZmaWXw4CM7OUcxCYmaWcg8DMLOUcBGZmKXdUvgvI5cQTT4xhw4bluwwzs8PGqlWr3oyIj3RnbEEGwbBhw6ivr893GWZmhw1Jr3V3rJeGzMxSzkFgZpZyDgIzs5QryNcIcnn//fdpbm5m586d+S7liFFcXExJSQl9+/bNdylmlkeHTRA0NzczYMAAhg0bhqR8l3PYiwi2bt1Kc3Mzw4cPz3c5ZpZHh83S0M6dOxk0aJBDoIdIYtCgQb7CMrNkQSBpoqT1kpokzTpAv7Ml7Zb0b10dm7COgxlu7Xg+zQwSBIGkIuAe4DKgFJgiqbSDfnOAx7o61szM8ifJawSVQFNEbACQtAioBta163cD8BBwdjfGWg4NDQ1s2rSJyy+/vEvjqqqqmDt3LhUVFb1UmdnBq62tZeHChfkuw0gWBEOA17P2m4FzsjtIGgJcDVzEvkHQ6disc9QANQCnnHJKgrKOfA0NDdTX13c5CMwOBwsXLmT5inr6DfabFfItSRDkWkhu/7VmdwNfj4jd7dadk4xtOxhRC9QCVFRUFOTXpv3yl79k3rx57Nq1i3POOYcvfOELzJgxg5UrV7J7924qKytZvHgxb775JrNnz2bQoEGsX7+e8ePH85Of/IQ+ffrw+OOP881vfpP33nuP0047jXvvvZdjjz2Wuro6brzxRv75z39y9NFH88QTTzB79mxaWlp47rnnuOWWW7jiiiu44YYbWL16Na2trXzrW9+iurqalpYWpk+fzrp16zjjjDNoaWnJ91SZJdJv8HA+OvWufJdxRHhtzhXdHpskCJqBoVn7JcCmdn0qgEWZEDgRuFxSa8KxXXbTTTfR0NBwsKfZR3l5OXfffXeH7S+99BKLFy/m+eefp2/fvnzpS19i/fr1XHXVVdx22220tLTw2c9+llGjRvHMM8+wcuVK1q1bx8c+9jEmTpzIww8/TFVVFXfeeSdPPvkk/fv3Z86cOXz/+99n1qxZTJo0icWLF3P22WfzzjvvcMwxx3D77bdTX1/Pj3/8YwBuvfVWLrroIhYsWMDbb79NZWUll1xyCT/96U855phjaGxspLGxkTFjxvTo3JjZkS1JENQBIyQNB/4KTAamZneIiL3XdpLuAx6NiEckHdXZ2MPF73//e1atWsXZZ7etfLW0tDB48GBmz57N2WefTXFxMfPmzdvbv7KyklNPPRWAKVOm8Nxzz1FcXMy6desYN24cALt27eLcc89l/fr1nHzyyXvPPXDgwJw1PP744yxdupS5c+cCbW+p/ctf/sLy5cv5yle+AkBZWRllZWW9MwlmdkTqNAgiolXS9bS9G6gIWBARayXNzLTP7+rYgy36QP9z7y0RwbXXXst3v/vdfY7/7W9/Y8eOHbz//vvs3LmT/v37A/u/NVMSEcGll17Kr371q33aGhsbE72VMyJ46KGH+PjHP75fm98KambdlehzBBGxLCJOj4jTIuL/Zo7NzxUCEfH5iHjwQGMPRxdffDEPPvggmzdvBuCtt97itddeo6amhjvuuINp06bx9a9/fW//lStX8uqrr/LBBx+wePFizj//fMaOHcvzzz9PU1MTAO+++y6vvPIKI0eOZNOmTdTV1QGwfft2WltbGTBgANu3b997zk9+8pP86Ec/IqLtJZQ//vGPAIwfP54HHngAgDVr1tDY2Nj7E2JmR4zD5hYT+VZaWsqdd97JhAkT+OCDD+jbty/V1dUcddRRTJ06ld27d3Peeefx1FNP0adPH84991xmzZrF6tWrGT9+PFdffTV9+vThvvvuY8qUKbz33nsA3HnnnZx++uksXryYG264gZaWFj70oQ/x5JNPcuGFF3LXXXdRXl7OLbfcwje+8Q1uuukmysrKiAiGDRvGo48+yhe/+EWmT59OWVkZ5eXlVFZW5nm2zOxwoj3/uywkFRUV0f6LaV566SXOOOOMPFXUNc888wxz587l0UcfzXcpnTqc5tWOLFVVVazYsNXvGuohr825YlVEdOvDQ4fNvYbMzKx3eGmoF1RVVVFVVZXvMszMEvEVgZlZyjkIzMxSzkFgZpZyDgIzs5Q7bF8sHjbrNz16vo13ferA7Rs3csUVV7BmzZp9js+ePZvx48dzySWX5Bz3yCOPcPrpp1Na6q9hMLPC5CuCg3T77bd3GALQFgTr1vXM1y+0trb2yHnMzLI5CLpg9+7dzJgxgzPPPJMJEybQ0tLC5z//eR58sO2OGrNmzaK0tJSysjK+9rWv8cILL7B06VJuvvlmysvL+fOf/0xDQwNjx46lrKyMq6++mn/84x8A1NXVUVZWxrnnnsvNN9/MqFGjALjvvvv4zGc+w5VXXsmECRPYsWMHF198MWPGjGH06NH8+te/BtquWEaOHMl1113HqFGjmDZtGk8++STjxo1jxIgRrFy5Mj+TZmYFz0HQBX/605/48pe/zNq1azn++ON56KGH9ra99dZbLFmyhLVr19LY2Mhtt93Geeedx1VXXcX3vvc9GhoaOO2007jmmmuYM2cOjY2NjB49mm9/+9sATJ8+nfnz5/Piiy9SVFS0z+O++OKL3H///Tz11FMUFxezZMkS/vCHP/D000/z1a9+de+9h5qamrjxxhtpbGzk5ZdfZuHChTz33HPMnTuX73znO4duoszssOIg6ILhw4dTXl4OwCc+8Qk2bty4t23gwIEUFxdz3XXX8fDDD3PMMcfsN37btm28/fbbXHDBBQBce+21LF++nLfffpvt27dz3nnnATB16r536r700kv58Ic/DLTdgfTWW2+lrKyMSy65hL/+9a/8/e9/31vf6NGj6dOnD2eeeSYXX3wxkhg9evQ+tZqZZXMQdMHRRx+9d7uoqGifNfujjjqKlStX8ulPf5pHHnmEiRMnJj5vZ/d72nNra4AHHniALVu2sGrVKhoaGjjppJPYuXPnfvX16dNn736fPn38+oKZdchB0EN27NjBtm3buPzyy7n77rv3foNa9q2kjzvuOE444QSeffZZAH7xi19wwQUXcMIJJzBgwABWrFgBwKJFizp8nG3btjF48GD69u3L008/zWuvvda7T8zMjniH7dtHO3u756G2fft2qqur2blzJxHBD37wAwAmT57MjBkzmDdvHg8++CD3338/M2fO5N133+XUU0/l3nvvBeDnP/85M2bMoH///lRVVXHcccflfJxp06Zx5ZVXUlFRQXl5OSNHjjxkz9HMjky+DXWB2LFjB8ceeywAd911F2+88QY//OEPe/1xj/R5tcLl21D3rF6/DbWkiZLWS2qSNCtHe7WkRkkNkuolnZ/VtlHS6j1t3SkyDX7zm99QXl7OqFGjePbZZ7ntttvyXZKZpUSnS0OSioB7gEuBZqBO0tKIyP6U1O+BpRERksqA/wKy1ywujIg3e7DuI86kSZOYNGlSvsswsxRKckVQCTRFxIaI2AUsAqqzO0TEjvj/a0z9gV5ZbyrEZazDmefTzCBZEAwBXs/ab84c24ekqyW9DPwG+EJWUwCPS1olqaajB5FUk1lWqt+yZct+7cXFxWzdutW/vHpIRLB161aKi4vzXYqZ5VmSdw0px7H9fhtHxBJgiaTxwB3AnhvwjIuITZIGA09IejkilucYXwvUQtuLxe3bS0pKaG5uJldIWPcUFxdTUlKS7zLMLM+SBEEzMDRrvwTY1FHniFgu6TRJJ0bEmxGxKXN8s6QltC017RcEnenbty/Dhw/v6jAzM+tEkqWhOmCEpOGS+gGTgaXZHST9iyRltscA/YCtkvpLGpA53h+YAOx7H2czM8urTq8IIqJV0vXAY0ARsCAi1kqamWmfD3wauEbS+0ALMCnzDqKTaFsu2vNYCyPid730XMzMrBsSfbI4IpYBy9odm5+1PQeYk2PcBuCsg6zRzMx6ke81ZGaWcg4CM7OUcxCYmaWcg8DMLOUcBGZmKecgMDNLOQeBmVnKOQjMzFLOQWBmlnIOAjOzlHMQmJmlnIPAzCzlHARmZinnIDAzSzkHgZlZyjkIzMxSLlEQSJooab2kJkmzcrRXS2qU1CCpXtL5SceamVl+dRoEkoqAe4DLgFJgiqTSdt1+D5wVEeXAF4CfdWGsmZnlUZIrgkqgKSI2RMQuYBFQnd0hInZERGR2+wORdKyZmeVXku8sHgK8nrXfDJzTvpOkq4HvAoOBT3VlrPWO2tpaFi5cmO8yjihTp06lpqYm32WY9agkVwTKcSz2OxCxJCJGAv8K3NGVsQCSajKvL9Rv2bIlQVnWmYULF9LQ0JDvMo4YDQ0NDlY7IiW5ImgGhmbtlwCbOuocEcslnSbpxK6MjYhaoBagoqIiZ1hY1+0cOJSNY2/OdxlHhJ0b/F4HOzIluSKoA0ZIGi6pHzAZWJrdQdK/SFJmewzQD9iaZKyZmeVXp1cEEdEq6XrgMaAIWBARayXNzLTPBz4NXCPpfaAFmJR58Tjn2F56LmZm1g1JloaIiGXAsnbH5mdtzwHmJB1rZmaFw58sNjNLOQeBmVnKOQjMzFLOQWBmlnIOAjOzlHMQmJmlnIPAzCzlHARmZinnIDAzSzkHgZlZyjkIzMxSzkFgZpZyDgIzs5RzEJiZpZyDwMws5RwEZmYp5yAwM0u5REEgaaKk9ZKaJO33Dd6SpklqzPy8IOmsrLaNklZLapBU35PFm5nZwev0qyolFQH3AJcCzUCdpKURsS6r26vABRHxD0mXAbXAOVntF0bEmz1Yt5mZ9ZAkVwSVQFNEbIiIXcAioDq7Q0S8EBH/yOyuAEp6tkwzM+stSYJgCPB61n5z5lhH/h34bdZ+AI9LWiWppqNBkmok1Uuq37JlS4KyzMysJ3S6NAQox7HI2VG6kLYgOD/r8LiI2CRpMPCEpJcjYvl+J4yopW1JiYqKipznNzOznpfkiqAZGJq1XwJsat9JUhnwM6A6IrbuOR4RmzJ/bgaW0LbUZGZmBSJJENQBIyQNl9QPmAwsze4g6RTgYeBzEfFK1vH+kgbs2QYmAGt6qngzMzt4nS4NRUSrpOuBx4AiYEFErJU0M9M+H5gNDAJ+IgmgNSIqgJOAJZljRwELI+J3vfJMzMysW5K8RkBELAOWtTs2P2v7OuC6HOM2AGe1P25mZoXDnyw2M0s5B4GZWco5CMzMUs5BYGaWcg4CM7OUcxCYmaWcg8DMLOUcBGZmKecgMDNLOQeBmVnKOQjMzFLOQWBmlnIOAjOzlHMQmJmlnIPAzCzlHARmZinnIDAzS7lEQSBpoqT1kpokzcrRPk1SY+bnBUlnJR1rZmb51WkQSCoC7gEuA0qBKZJK23V7FbggIsqAO4DaLow1M7M8SvKdxZVAU+b7h5G0CKgG1u3pEBEvZPVfAZQkHWt2OFm+op7iU0bnu4wjwq7Nr9Jv8PB8l2EkWxoaAryetd+cOdaRfwd+29Wxkmok1Uuq37JlS4KyzA6t/qVV/sXVg/oNHk7/0qp8l2EkuyJQjmORs6N0IW1BcH5Xx0ZELZklpYqKipx9zPJpQPlEBpRPzHcZZj0uSRA0A0Oz9kuATe07SSoDfgZcFhFbuzLWzMzyJ8nSUB0wQtJwSf2AycDS7A6STgEeBj4XEa90ZayZmeVXp1cEEdEq6XrgMaAIWBARayXNzLTPB2YDg4CfSAJojYiKjsb20nMxM7NuSLI0REQsA5a1OzY/a/s64LqkY83MrHD4k8VmZinnIDAzSzkHgZlZyjkIzMxSzkFgZpZyDgIzs5RzEJiZpZyDwMws5RwEZmYp5yAwM0s5B4GZWco5CMzMUs5BYGaWcg4CM7OUcxCYmaWcg8DMLOUcBGZmKZcoCCRNlLReUpOkWTnaR0p6UdJ7kr7Wrm2jpNWSGiTV91ThZmbWMzr9qkpJRcA9wKVAM1AnaWlErMvq9hbwFeBfOzjNhRHx5kHWamZmvSDJFUEl0BQRGyJiF7AIqM7uEBGbI6IOeL8XajQzs16UJAiGAK9n7TdnjiUVwOOSVkmq6aiTpBpJ9ZLqt2zZ0oXTm5nZwUgSBMpxLLrwGOMiYgxwGfBlSeNzdYqI2oioiIiKj3zkI104vZmZHYwkQdAMDM3aLwE2JX2AiNiU+XMzsIS2pSYzMysQSYKgDhghabikfsBkYGmSk0vqL2nAnm1gArCmu8WamVnP6/RdQxHRKul64DGgCFgQEWslzcy0z5f0UaAeGAh8IOkmoBQ4EVgiac9jLYyI3/XKMzEzs27pNAgAImIZsKzdsflZ23+jbcmovXeAsw6mQDMz613+ZLGZWco5CMzMUs5BYGaWcg4CM7OUcxCYmaWcg8DMLOUcBGZmKecgMDNLOQeBmVnKOQjMzFLOQWBmlnIOAjOzlHMQmJmlnIPAzCzlHARmZinnIDAzSzkHgZlZyiUKAkkTJa2X1CRpVo72kZJelPSepK91ZayZmeVXp0EgqQi4B7iMtu8hniKptF23t4CvAHO7MdbMzPIoyXcWVwJNEbEBQNIioBpYt6dDRGwGNkv6VFfH5rJ+/XqqqqqSPgfrQENDAwwcmu8yzKzAJQmCIcDrWfvNwDkJz594rKQaoAaAor6s2LA14UNYhwYOpX9pVb6rMLMClyQIlONYJDx/4rERUQvUAhx98oj46NS7Ej6EmZkdjCQvFjcD2esLJcCmhOc/mLFmZnYIJAmCOmCEpOGS+gGTgaUJz38wY83M7BDodGkoIlolXQ88BhQBCyJiraSZmfb5kj4K1AMDgQ8k3QSURsQ7ucb20nMxM7NuSPIaARGxDFjW7tj8rO2/0bbsk2ismZkVDn+y2Mws5RwEZmYp5yAwM0s5B4GZWco5CMzMUs5BYGaWcg4CM7OUcxCYmaWcg8DMLOUcBGZmKecgMDNLOQeBmVnKOQjMzFLOQWBmlnIOAjOzlHMQmJmlnIPAzCzlEgWBpImS1ktqkjQrR7skzcu0N0oak9W2UdJqSQ2S6nuyeDMzO3idflWlpCLgHuBSoBmok7Q0ItZldbsMGJH5OQf4j8yfe1wYEW/2WNVmZtZjklwRVAJNEbEhInYBi4Dqdn2qgf+MNiuA4yWd3MO1mplZL0gSBEOA17P2mzPHkvYJ4HFJqyTVdPQgkmok1Uuq3/3utgRlmZlZT+h0aQhQjmPRhT7jImKTpMHAE5Jejojl+3WOqAVqAY4+eUT785uZWS9JckXQDAzN2i8BNiXtExF7/twMLKFtqcnMzApEkiCoA0ZIGi6pHzAZWNquz1Lgmsy7h8YC2yLiDUn9JQ0AkNQfmACs6cH6zczsIHW6NBQRrZKuBx4DioAFEbFW0sxM+3xgGXA50AS8C0zPDD8JWCJpz2MtjIjf9fizMDOzbkvyGgERsYy2X/bZx+ZnbQfw5RzjNgBnHWSNZmbWi/zJYjOzlHMQmJmlnIPAzCzlHARmZinnIDAzSzkHgZlZyjkIzMxSzkFgZpZyDgIzs5RzEJiZpZyDwMws5RwEZmYp5yAwM0s5B4GZWco5CMzMUs5BYGaWcg4CM7OUSxQEkiZKWi+pSdKsHO2SNC/T3ihpTNKxZmaWX50GgaQi4B7gMqAUmCKptF23y4ARmZ8a4D+6MNbMzPIoyRVBJdAUERsiYhewCKhu16ca+M9oswI4XtLJCceamVkeJfny+iHA61n7zcA5CfoMSTgWAEk1tF1NALz32pwr1iSoLZ9OBN7MdxEJuM6e5Tp7luvsOR/v7sAkQaAcxyJhnyRj2w5G1AK1AJLqI6IiQW15czjUCK6zp7nOnuU6e46k+u6OTRIEzcDQrP0SYFPCPv0SjDUzszxK8hpBHTBC0nBJ/YDJwNJ2fZYC12TePTQW2BYRbyQca2ZmedTpFUFEtEq6HngMKAIWRMRaSTMz7fOBZcDlQBPwLjD9QGMT1FXbnSdziB0ONYLr7Gmus2e5zp7T7RoVkXPJ3szMUsKfLDYzSzkHgZlZyuUtCBLctqJK0jZJDZmf2Xmqc4GkzZJyfq7hQLfXKKAaC2Uuh0p6WtJLktZKujFHn0KYzyR15n1OJRVLWinpfzJ1fjtHn7zOZ8Ia8z6XWbUUSfqjpEdztOX932ZWLQeqs+vzGRGH/Ie2F47/DJxK21tM/wcobdenCng0H/W1q2M8MAZY00H75cBvafvMxFjgvwuwxkKZy5OBMZntAcArOf7eC2E+k9SZ9znNzNGxme2+wH8DYwtpPhPWmPe5zKrlfwMLc9WT77nsQp1dns98XREcNreeiIjlwFsH6NLR7TUOmQQ1FoSIeCMi/pDZ3g68RNunz7MVwnwmqTPvMnO0I7PbN/PT/t0feZ3PhDUWBEklwKeAn3XQJe//NiFRnV2WryDo6JYU7Z2buaT8raQzD01pXZb0ueRbQc2lpGHA/6Ltf4jZCmo+D1AnFMCcZpYIGoDNwBMRUXDzmaBGKIC5BO4G/g/wQQfteZ/LjLs5cJ3QxfnMVxAkufXEH4CPRcRZwI+AR3q7qG5KfBuNPCqouZR0LPAQcFNEvNO+OceQvMxnJ3UWxJxGxO6IKKftU/uVkka165L3+UxQY97nUtIVwOaIWHWgbjmOHdK5TFhnl+czX0HQ6W0rIuKdPZeUEbEM6CvpxENXYmJJbsGRV4U0l5L60vbL9YGIeDhHl4KYz87qLKQ5zdTwNvAMMLFdU0HMJ3RcY4HM5TjgKkkbaVuqvkjSL9v1KYS57LTO7sxnvoKg01tPSPqoJGW2K2mrdeshr7RzHd1eo2AUylxmavg58FJEfL+DbnmfzyR1FsKcSvqIpOMz2x8CLgFebtctr/OZpMZCmMuIuCUiSiJiGG2/j56KiM+265b3f5tJ6uzOfCa56VyPi2S3rfg34IuSWoEWYHJkXhI/lCT9irZX4U+U1Ax8k7YXvPbUmfP2GgVWY0HMJW3/m/kcsDqzZgxwK3BKVq15n0+S1VkIc3oycL/avgCqD/BfEfGoEtz+pcBqLIS5zKnA5rJDBzufvsWEmVnK+ZPFZmYp5yAwM0s5B4GZWco5CMzMUs5BYGaWcg4CM7OUcxCYmaXc/wMBrac8GGFOJgAAAABJRU5ErkJggg==\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAWsElEQVR4nO3de5hWZb3/8fd3hoMOKppYmSCQoYSHklBE1I3u3CJbRdtYGD/bbW2bpqblMfNnhnncam4FRUy2uX952m4pTBRNLc0DgYQKKokkMlGZhiBnBu7fHzPiMA7OM7pmnkfu9+u65mId7rWe77Oumw+LdbgnUkpIkvJSVe4CJEntz/CXpAwZ/pKUIcNfkjJk+EtShjqU64O7deuWevXqVa6Pl6SPpGeeeeaNlNJ2H3Y/ZQv/Xr16MX369HJ9vCR9JEXE/CL242UfScqQ4S9JGTL8JSlDhr8kZcjwl6QMGf6SlKGSwj8ihkbEnIiYGxHnNrN+SEQsjoiZDT8XFF+qJKkoLT7nHxHVwFjgYKAWmBYRk1JKLzRp+nhK6bA2qFGSVLBSXvLaG5ibUpoHEBF3AMOBpuHfKnPmzGHIkCEfZheSpA+olPDfAVjQaL4WGNhMu0ER8SywEDgzpTS7aYOIOAE4AYDqjjw9781WFyxJ+vBKCf9oZlnTX/81A+iZUloaEcOAnwN93rNRSuOB8QCdt++TPvnVy1pXrSRlbv7lxVxdL+WGby3Qo9F8d+rP7tdLKS1JKS1tmJ4MdIyIboVUKEkqXCnhPw3oExG9I6ITMBKY1LhBRHwyIqJheu+G/XpNR5IqVIuXfVJKdRFxCjAFqAYmpJRmR8SJDevHASOAkyKiDlgBjEz+ZnhJqlglDenccClncpNl4xpNjwHGFFuaJKmt+IavJGXI8JekDBn+kpQhw1+SMmT4S1KGDH9JypDhL0kZMvwlKUOGvyRlyPCXpAwZ/pKUIcNfkjJk+EtShgx/ScqQ4S9JGTL8JSlDhr8kZcjwl6QMGf6SlCHDX5IyZPhLUoYMf0nKkOEvSRky/CUpQ4a/JGXI8JekDBn+kpQhw1+SMmT4S1KGDH9JypDhL0kZMvwlKUMlhX9EDI2IORExNyLOfZ92e0XE2ogYUVyJkqSitRj+EVENjAUOBfoBx0REv420uxyYUnSRkqRilXLmvzcwN6U0L6W0GrgDGN5Mu1OB/wVeL7A+SVIbKCX8dwAWNJqvbVi2XkTsABwFjHu/HUXECRExPSKmr12+uLW1SpIKUkr4RzPLUpP5a4BzUkpr329HKaXxKaUBKaUB1TVdSyxRklS0DiW0qQV6NJrvDixs0mYAcEdEAHQDhkVEXUrp50UUKUkqVinhPw3oExG9gT8BI4GvNm6QUur9znRE3AL80uCXpMrVYvinlOoi4hTqn+KpBiaklGZHxIkN69/3Or8kqfKUcuZPSmkyMLnJsmZDP6X09Q9fliSpLfmGryRlyPCXpAwZ/pKUIcNfkjJk+EtShgx/ScqQ4S9JGTL8JSlDhr8kZcjwl6QMGf6SlCHDX5IyZPhLUoYMf0nKkOEvSRky/CUpQ4a/JGXI8JekDBn+kpQhw1+SMmT4S1KGDH9JypDhL0kZMvwlKUOGvyRlyPCXpAwZ/pKUIcNfkjJk+EtShgx/ScqQ4S9JGTL8JSlDJYV/RAyNiDkRMTcizm1m/fCIeC4iZkbE9IjYr/hSJUlF6dBSg4ioBsYCBwO1wLSImJRSeqFRs4eBSSmlFBF7AHcBfduiYEnSh1fKmf/ewNyU0ryU0mrgDmB44wYppaUppdQw2wVISJIqVinhvwOwoNF8bcOyDUTEURHxEnAfcFwx5UmS2kIp4R/NLHvPmX1KaWJKqS9wJHBRszuKOKHhnsD0tcsXt6pQSVJxSgn/WqBHo/nuwMKNNU4pPQbsFBHdmlk3PqU0IKU0oLqma6uLlSQVo5Twnwb0iYjeEdEJGAlMatwgIj4TEdEw3R/oBLxZdLGSpGK0+LRPSqkuIk4BpgDVwISU0uyIOLFh/TjgX4CvRcQaYAXwlUY3gCVJFabF8AdIKU0GJjdZNq7R9OXA5cWWJklqK77hK0kZMvwlKUOGvyRlyPCXpAwZ/pKUIcNfkjJk+EtShgx/ScqQ4S9JGTL8JSlDhr8kZcjwl6QMGf6SlCHDX5IyZPhLUoYMf0nKkOEvSRky/CUpQ4a/JGXI8JekDBn+kpQhw1+SMmT4S1KGDH9JypDhL0kZMvwlKUOGvyRlyPCXpAwZ/pKUIcNfkjJk+EtShgx/ScqQ4S9JGSop/CNiaETMiYi5EXFuM+tHRcRzDT9PRsTnii9VklSUFsM/IqqBscChQD/gmIjo16TZH4F/SCntAVwEjC+6UElScUo5898bmJtSmpdSWg3cAQxv3CCl9GRKaVHD7NNA92LLlCQVqZTw3wFY0Gi+tmHZxhwP3N/ciog4ISKmR8T0tcsXl16lJKlQHUpoE80sS802jDiQ+vDfr7n1KaXxNFwS6rx9n2b3IUlqe6WEfy3Qo9F8d2Bh00YRsQfwE+DQlNKbxZQnSWoLpVz2mQb0iYjeEdEJGAlMatwgInYE7gGOTSn9ofgyJUlFavHMP6VUFxGnAFOAamBCSml2RJzYsH4ccAGwLXB9RADUpZQGtF3ZkqQPo5TLPqSUJgOTmywb12j6G8A3ii1NktRWfMNXkjJk+EtShgx/ScqQ4S9JGTL8JSlDhr8kZcjwl6QMGf6SlCHDX5IyZPhLUoYMf0nKkOEvSRky/CUpQ4a/JGXI8JekDBn+kpQhw1+SMmT4S1KGDH9JypDhL0kZMvwlKUOGvyRlyPCXpAwZ/pKUIcNfkjJk+EtShgx/ScqQ4S9JGTL8JSlDhr8kZcjwl6QMGf6SlKGSwj8ihkbEnIiYGxHnNrO+b0Q8FRGrIuLM4suUJBWpQ0sNIqIaGAscDNQC0yJiUkrphUbN/g58GziyLYqUJBWrlDP/vYG5KaV5KaXVwB3A8MYNUkqvp5SmAWvaoEZJUsFKCf8dgAWN5msblrVaRJwQEdMjYvra5Ys/yC4kSQUoJfyjmWXpg3xYSml8SmlASmlAdU3XD7ILSVIBSgn/WqBHo/nuwMK2KUeS1B5KCf9pQJ+I6B0RnYCRwKS2LUuS1JZafNonpVQXEacAU4BqYEJKaXZEnNiwflxEfBKYDmwFrIuI04F+KaUlrSlmq85VnDpwG3pu3ZFo9mqTWiuRmP/WGq6buoglq9aVuxxJFaLF8AdIKU0GJjdZNq7R9F+ovxz0oZw6cBv67/QpOtRsSYThX4SUEttuu4RTgYsfe7Pc5UiqEBX1hm/PrTsa/AWLCDrUbEXPrTuWuxRJFaSiwj8Ig78NRISX0SRtoKLCX5LUPgz/CvDS7Od5/JEHW73d8Ucfxuxnf98GFUna1Bn+FWDO7Od5/JGHyl2GpIyU9LRPOVxx4feYM/v5Qve5y667c/aFl7bY7pf33MltE8ZTt2Y1u+35BY78yv9h9Nmn8bN7H2bturWMOvyLXDH2Zt5a9Heuv/ISum7zMebPm0v/gYP4/sVXUVVVxZO/eYQbrr6M1atX0aNnb0ZfNYaaLlswa+YMrrjwXFYsX07HTp258bZ7uP6qS1i1ciUzpz3NcSd/hwO+eAiX/d9zmPvSC9StreOk75zLgYcMY+WKFVxwxsnMe3kOvT+zCytXriz0+EjKR8WGf7nMe3kOU+6dyE8nPkDHjh25+LwzmP/KXIYcfChj/uNiVq1cwWFHHU2fvv2Y9tRvmfXsDCY+/DTbd+/Bt44dwcP338uAQftx07VXcuPtE6mp6cKE66/h1puu5/hvnc7ZJx/HFWMnsNvn+7P07SVstnkN3zrjPGY/93vO+9F/AHDtZaPZe/D+jL5qDEsWL2bU4f/IwP3/gbv/3y1stnkNdz/0BH94cRYjDx1S3oMl6SOrYsO/lDP0tjD1id/w4nPPMuqwgwBYuXIlH+u2Hd88/Wy+ethBdOq8GeeMvnx9+90+15/uPXsBMPSIf+H3056mU+fOzHt5Dl8/aigAa9asYY/+e/HqKy+z3cc/wW6f7w/AFltu1WwNTz32KL9+6H5uvXEMAKtXreQvf6plxtQnOea4bwKw82d3o89nd22TYyBp01ex4V8uKcHhR4/ktHN/sMHyN17/K8uXLaOuro5Vq1ZSU9MF4D2PpgZBSol99h/C5WNv3mDdH16cBSU8yppIXD3+Vnrt1Oc963wUVlIRvOHbxMDBB/Cr+ybx5ht/A2DxokUsrH2N0eeczslnnsewI0dwzSUXrm8/a+YMal+bz7p165hy70T23Hsf9ui/FzOnT+W1P84DYMWK5bw6by69d9qZv/31L8yaOQOAZUvfpq6uji5bbMHyZUvX73PfAw7itv8aT0r1g6e+OOs5APoP3JfJE/8HgJdfeoGXX5zd5sdD0qbJM/8mdtq5Lyef9X1OGvUl1q1bR4eOHRnyT8Oo7tCBYUcdzdq1a/nakYcw9YnHqKqqYo8v7MV/XvpD5s55gf4DB3HQ0MOoqqpi9NXXc+4p32D16lUAnHLW9+n16c9wxdgJXHbBOaxauYLOm23O+Nsnsteg/Zkw9hq+fMj+HHfydzjhtLO44offY8TBg0kp8akeOzLmljv58rHHccEZJzPi4MHssuvu6y8fSVJrxTtnl+2t8/Z90vb/es0Gy246Yns+seOny1LPBzHtqd/y0xuvY8wtd5a7lBb99bV5/PukP5e7DEkf0vzLD3smpTTgw+7Hyz6SlCEv+3wIew3aj70G7VfuMiSp1Tzzl6QMGf6SlCHDX5IyZPhLUoYq+obvEWOeKHR/k04Z3GKbPy14jVO//hXuefipDZaPvfISvjBwX/bZf0iz2z3ywH30/PRO7LRz3yJKlaQ25Zl/iU4+87yNBj/Ao1PuY97Lcwr5rLq6ukL2I0kbU9Fn/uWybt06fnj2acx85nd8/BPb8583/4yLzzuDA754CAf/83CuufRCfvPQA1RXVzPogIP4x0MP49cP3c/0qU9w07VXctWNt7Js2VJ+9L3vsnLFcrr37M3oK8ew1dZbM2vmDC4861Q2r+nCnnsN5LeP/op7Hn6KX9x1G48/8iCrVq1kxfLlXDvhNk47fhRLFr9F3Zo1nHLW+Rx4yDD+tOA1vnXsCPbcax+emzGNXfrtxvAvj+KGqy/l72+8wSXXjmf3Pb9Q7kMoqcJ55t+M1/74Cl/5128w8eGn2KprV351/6T16xYvWsQjD9zHPQ8/xd0PPcG/f/tMPj9gIEMOPpTvfn80d015nB69enP+6Sdy+vcu5O6HnqBP336Mu6Z+JNALzjiZ8y+9mv/+xYNUVVdv8LnPPjONH119Az+5cxKdOm/Gj2/6b+68/zf85K57ueqi89eP9bPg1XmMOu6b3P3QE/zxlZeZ/PO7ueWeB/ju+Rdx85ir2+9ASfrIMvybsUOPnvTddXcAPrv751i4YMH6dV223JLOnTtz4Vnf5lf338vmm2/+nu3fXrKYt5csZsCg+nsMR4w4hmemPsmSxYtZvmwpnx8wEIBhR47YYLt99h9C1222ASClxLWXX8SIgwfzzWOO5PW//Jk3//b6+vr6fHZXqqqq2Gnnvgzc7wAigj59+7Gw9rXiD4ikTY7h34yOnTqtn66uqmbt2nevwXfo0IGf3fswXxx2OI9OuY+Tjh3R3C6a18I4SpvX1Kyfnjzxf1j05pvcPvnX3DXlcbbdbjtWrVr1nvqqoopOnToDEFVV3i+QVBLDv5WWL1vK228vYf+D/omzf3Dp+l81WbPFFixbWj8s85ZbdWWrrlszY+qTQP2vhRwwcDBbbb01NV224LkZ0wB44Bf3bPRzlr69hI9160bHjh353ZOPs7B2wUbbSlJrVfQN31IezWxvy5Yu5bTjR7F61UpSSpz1g0sAGHrElxh9zunc9l83ctW4n3LRj29494bvjr0YfdVYAC688jpGn30am9d0YcCgwWy5VfO/zWvYUUfz7X87hmOGHcguu+5O78/s3G7fUdKmzyGd29nyZUup6bIFADeP/TFvvP5XzvnhZW3+uQ7pLG0aihrSuaLP/DdFjz38IBPG/pi6ujo+1b0Ho6++vtwlScqQ4d/Ohh7xJYYe8aVylyEpcxV1wzeRKNdlqE1ZSomEx1XSuyoq/Oe/tYa65Uv8B6BAKSXqli9h/ltryl2KpApSUZd9rpu6iFOBnlu/QRDlLmeTkEjMf2sN101dVO5SJFWQigr/JavWcfFjb5a7DEna5JV02ScihkbEnIiYGxHnNrM+IuLahvXPRUT/4kuVJBWlxfCPiGpgLHAo0A84JiL6NWl2KNCn4ecE4IaC65QkFaiUM/+9gbkppXkppdXAHcDwJm2GA7emek8DW0fE9gXXKkkqSCnX/HcAGg8sUwsMLKHNDsAGr5RGxAnU/88AYNX8yw+b1apqy6Mb8Ea5iyiBdRbro1DnR6FGsM6i7VLETkoJ/+Yeu2n6LGYpbUgpjQfGA0TE9CJeUW5r1lks6yzOR6FGsM6iRcT0IvZTymWfWqBHo/nuwMIP0EaSVCFKCf9pQJ+I6B0RnYCRwKQmbSYBX2t46mcfYHFKyVHEJKlCtXjZJ6VUFxGnAFOAamBCSml2RJzYsH4cMBkYBswFlgP/VsJnj//AVbcv6yyWdRbno1AjWGfRCqmzbEM6S5LKp6LG9pEktQ/DX5Iy1CbhX8JwEEMiYnFEzGz4uaDUbduxxrMa1TcrItZGxMca1r0aEc83rCvksav3qXNCRLweEc2+E/F+Q2u017Essc5RDfU9FxFPRsTnGq2rpONZCX2zpRorpW/2iIhHI+LFiJgdEac106bs/bPEOsveP0uss7j+mVIq9If6m8KvAJ8GOgHPAv2atBkC/PKDbNteNTZpfzjwSKP5V4FuRde1kc8+AOgPzNrI+mHA/dS/a7EPMLU9j2Ur6twX2KZh+tB36qzA41nWvllKjRXUN7cH+jdMbwn8oZm/62XvnyXWWfb+WWKdhfXPtjjzL2U4iLbYti1rPAa4vQ3qaFFK6THg7+/TZGNDa7TXsSypzpTSkymld8aVfpr6d0HaXQnHc2Pa7Xi2ssZy9s0/p5RmNEy/DbxI/Zv9jZW9f5ZSZyX0zxKP58a0+ni2RfhvbKiHpgZFxLMRcX9E7NrKbdurRiKiBhgK/G+jxQl4MCKeifohK8ppY9+lvY7lB3E89WeD76ik4wnl7Zslq6S+GRG9gD2BqU1WVVT/fJ86Gyt7/2yhzkL6Z1uM51/KUA8zgJ4ppaURMQz4OfUjgpY0TEQBWvM5hwNPpJQan4kNTiktjIiPAw9FxEsNZ2vlsLHv0l7HslUi4kDq/3Lt12hxJR3PcvfN1qiIvhkRW1D/D9DpKaUlTVc3s0lZ+mcLdb7Tpuz9s4U6C+ufbXHm3+JQDymlJSmlpQ3Tk4GOEdGtlG3bq8ZGRtLkv9UppYUNf74OTKT+v1zlsrHvUnFDbkTEHsBPgOEppfW/taeSjmcF9M3WKHvfjIiO1AfVz1JK9zTTpCL6Zwl1VkT/bKnOQvtnG9y06ADMA3rz7o2HXZu0+STvvmC2N/Aa9f9ytbhte9XY0K4r9ddeuzRa1gXYstH0k8DQomtsUkcvNn6D8p/Z8Iba71rzHduxzh2pfwN83ybLK+14lrVvllJjpfTNhuNyK3DN+7Qpe/8ssc6y988S6yysfxZ+2SeVNhzECOCkiKgDVgAjU/23aXbbMtUIcBTwYEppWaPNPwFMjAioP+C3pZQeKLrGd0TE7dTf4e8WEbXAD4COjepsdmiNjX3HMtZ5AbAtcH3DsatL9SMoVtrxLGvfLLFGqIC+CQwGjgWej4iZDcvOoz5IK6l/llJnJfTPUuosrH86vIMkZcg3fCUpQ4a/JGXI8JekDBn+kpQhw1+SMmT4S1KGDH9JytD/BxNLwX89BbUSAAAAAElFTkSuQmCC\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAX7ElEQVR4nO3df3BV5b3v8fc3IRBBQAulMoQj0YvFCGkuhQjoxbT8aOSgOUzb4VePlh5hsMXqTOs0Opa26q0wMq3FH4dmThFvK4UziDSD9GipcigqhwSbRn4Ym6JoRItCQZAgBr/3j2yZbdwxK8ne7Ow8n9dMZtaPZ6/9ffbDfFh59lor5u6IiEj3l5XuAkRE5OxQ4IuIBEKBLyISCAW+iEggFPgiIoHoka43HjhwoA8bNixdby8ikpF27tz5jrt/tiOvTVvgDxs2jOrq6nS9vYhIRjKz/R19raZ0REQCocAXEQmEAl9EJBBpm8MXke7hgw8+oKGhgZMnT6a7lG4lNzeXvLw8cnJyknZMBb6IdEpDQwN9+/Zl2LBhmFm6y+kW3J1Dhw7R0NBAfn5+0o6rKR0R6ZSTJ08yYMAAhX0SmRkDBgxI+m9NkQLfzErNrM7M6s2sPMH+W82sJvazy8xOm9lnklqpiHRZCvvkS8Vn2mbgm1k28CBwNVAAzDazgvg27n6vuxe5exFwG/Df7n446dWKiEiHRZnDLwbq3X0fgJmtAcqAPa20nw38tq2D1tXVUVJSErFM6WrmzJnDggUL0l2GSNrU1NRw4MABpk2b1q7XlZSUsGzZMsaMGZOiyloXJfCHAK/HrTcAlydqaGa9gVJgUSv7FwDNKZGdw/Z9h9pTq3QRpw6+AqDAl6DV1NRQXV3d7sBPpyiBn2giqbU/k3UN8Gxr0znuXgFUAPQaPNwvmLMkUpHStby1+hNf44ik3W9+8xuWL1/OqVOnuPzyy/nWt77F/Pnz2bFjB6dPn6a4uJi1a9fyzjvvsHjxYgYMGEBdXR0TJ07koYceIisri6eeeoof/ehHvP/++1x88cU8/PDDnHvuuVRVVXHzzTfz3nvv0atXL/7whz+wePFiGhsb2bZtG7fddhvTp0/npptu4sUXX6SpqYkf//jHlJWV0djYyLx589izZw+XXnopjY2NafuMogR+AzA0bj0PONBK21lEmM4Rke7plltuoaamJqnHLCoq4r777vvUNnv37mXt2rU8++yz5OTk8O1vf5u6ujquvfZa7rjjDhobG/nGN77ByJEj2bJlCzt27GDPnj1ceOGFlJaWsn79ekpKSrj77rvZvHkzffr0YenSpfzsZz+jvLycmTNnsnbtWsaOHcu7775L7969ufPOO6muruaBBx4A4Pbbb+fLX/4yK1eu5MiRIxQXFzN58mR++ctf0rt3b2pra6mtrWX06NFJ/XzaI0rgVwHDzSwfeIPmUJ/TspGZ9QeuAr6R1ApFRNrwxz/+kZ07dzJ27FgAGhsbGTRoEIsXL2bs2LHk5uayfPnyM+2Li4u56KKLAJg9ezbbtm0jNzeXPXv2cMUVVwBw6tQpxo8fT11dHYMHDz5z7H79+iWs4amnnqKyspJly5YBzZervvbaa2zdupXvfve7ABQWFlJYWJiaDyGCNgPf3ZvMbBHwJJANrHT33Wa2MLZ/RazpDOApd38vZdWKSJfW1pl4qrg7119/Pffcc8/Htr/11lscP36cDz74gJMnT9KnTx/gk5c8mhnuzpQpU/jtbz8+SVFbWxvpEkl357HHHuPzn//8J/Z1lctWI12H7+6b3P0Sd7/Y3f9vbNuKuLDH3Ve5+6xUFSoi0ppJkyaxbt06Dh48CMDhw4fZv38/CxYs4K677mLu3Ln84Ac/ONN+x44dvPLKK3z44YesXbuWK6+8knHjxvHss89SX18PwIkTJ3j55ZcZMWIEBw4coKqqCoBjx47R1NRE3759OXbs2JljfuUrX+H+++/Hvfkrzj//+c8ATJw4kUcffRSAXbt2UVtbm/oPpBV6tIKIZLyCggLuvvtupk6dyocffkhOTg5lZWX06NGDOXPmcPr0aSZMmMDTTz9NVlYW48ePp7y8nBdffJGJEycyY8YMsrKyWLVqFbNnz+b9998H4O677+aSSy5h7dq13HTTTTQ2NnLOOeewefNmvvSlL7FkyRKKioq47bbb+OEPf8gtt9xCYWEh7s6wYcPYuHEjN954I/PmzaOwsJCioiKKi4vT9jnZR/8bnW29Bg/3wdffl5b3ls55a3U54y4awJYtW9JdinQBe/fu5dJLL013GZFt2bKFZcuWsXHjxnSX0qZEn62Z7XT3Dl3Er2fpiIgEQlM6IhKUkpKSYO/y1xm+iEggFPgiIoFQ4IuIBEKBLyISCH1pKyJJNaz8iaQe79Ul/9x2m1dfZfr06ezatetj2xcvXszEiROZPHlywtdt2LCBSy65hIKCgoT7uxud4YtIt3XnnXe2GvbQHPh79rT2pz3ap6mpKSnHSSUFvoh0C6dPn2b+/PlcdtllTJ06lcbGRr75zW+ybt06AMrLyykoKKCwsJDvf//7PPfcc1RWVnLrrbdSVFTE3/72N2pqahg3bhyFhYXMmDGDf/zjHwBUVVVRWFjI+PHjufXWWxk5ciQAq1at4utf/zrXXHMNU6dO5fjx40yaNInRo0czatQofve73wHNv4GMGDGCG264gZEjRzJ37lw2b97MFVdcwfDhw9mxY8dZ+YwU+CLSLfz1r3/lO9/5Drt37+a8887jscceO7Pv8OHDPP744+zevZva2lruuOMOJkyYwLXXXsu9995LTU0NF198Mddddx1Lly6ltraWUaNG8ZOf/ASAefPmsWLFCp5//nmys7M/9r7PP/88jzzyCE8//TS5ubk8/vjjvPDCCzzzzDN873vfO/Nsnfr6em6++WZqa2t56aWXWL16Ndu2bWPZsmX89Kc/PSufkQJfRLqF/Px8ioqKAPjiF7/Iq6++emZfv379yM3N5YYbbmD9+vX07t37E68/evQoR44c4aqrrgLg+uuvZ+vWrRw5coRjx44xYcIEoPnPe8abMmUKn/nMZ4DmJ2befvvtFBYWMnnyZN544w3+/ve/n6lv1KhRZGVlcdlllzFp0iTMjFGjRn2s1lRS4ItIt9CrV68zy9nZ2R+bU+/Rowc7duzgq1/9Khs2bKC0tDTycdt63thHj1wGePTRR3n77bfZuXMnNTU1fO5zn+PkyZOfqC8rK+vMelZW1lmb/1fgi0i3d/z4cY4ePcq0adO47777zvxVrvhHHPfv35/zzz+fP/3pTwD8+te/5qqrruL888+nb9++bN++HYA1a9a0+j5Hjx5l0KBB5OTk8Mwzz7B///7UdqyddFmmiCRVlMsoz7Zjx45RVlbGyZMncXd+/vOfAzBr1izmz5/P8uXLWbduHY888ggLFy7kxIkTXHTRRTz88MMA/OpXv2L+/Pn06dOHkpIS+vfvn/B95s6dyzXXXMOYMWMoKipixIgRZ62PUejxyNJuejyyxMu0xyN3xPHjxzn33HMBWLJkCW+++Sa/+MUvUv6+yX48ss7wRUTa8MQTT3DPPffQ1NTEhRdeyKpVq9JdUoco8EVE2jBz5kxmzpyZ7jI6TV/aikinpWtquDtLxWcaKfDNrNTM6sys3szKW2lTYmY1ZrbbzP47uWWKSFeVm5vLoUOHFPpJ5O4cOnSI3NzcpB63zSkdM8sGHgSmAA1AlZlVuvueuDbnAQ8Bpe7+mpkNSmqVItJl5eXl0dDQwNtvv53uUrqV3Nxc8vLyknrMKHP4xUC9u+8DMLM1QBkQ/8ShOcB6d38NwN0PJrVKEemycnJyyM/PT3cZEkGUKZ0hwOtx6w2xbfEuAc43sy1mttPMrkt0IDNbYGbVZlZ9+sTRjlUsIiIdEuUM3xJsazlZ1wP4IjAJOAd43sy2u/vLH3uRewVQAc3X4be/XBER6agogd8ADI1bzwMOJGjzjru/B7xnZluBLwAvIyIiXUKUKZ0qYLiZ5ZtZT2AWUNmize+A/2NmPcysN3A5sDe5pYqISGe0eYbv7k1mtgh4EsgGVrr7bjNbGNu/wt33mtl/AbXAh8B/uPuu1o8qIiJnW6Q7bd19E7CpxbYVLdbvBe5NXmkiIpJMutNWRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQkQLfzErNrM7M6s2sPMH+EjM7amY1sZ/FyS9VREQ6o0dbDcwsG3gQmAI0AFVmVunue1o0/ZO7T09BjSIikgRRzvCLgXp33+fup4A1QFlqyxIRkWSLEvhDgNfj1hti21oab2Z/MbPfm9lliQ5kZgvMrNrMqk+fONqBckVEpKOiBL4l2OYt1l8ALnT3LwD3AxsSHcjdK9x9jLuPye7dv12FiohI50QJ/AZgaNx6HnAgvoG7v+vux2PLm4AcMxuYtCpFRKTTogR+FTDczPLNrCcwC6iMb2BmF5iZxZaLY8c9lOxiRUSk49q8Ssfdm8xsEfAkkA2sdPfdZrYwtn8F8DXgRjNrAhqBWe7ectpHRETSqM3AhzPTNJtabFsRt/wA8EBySxMRkWTSnbYiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBEKBLyISCAW+iEggIgW+mZWaWZ2Z1ZtZ+ae0G2tmp83sa8krUUREkqHNwDezbOBB4GqgAJhtZgWttFsKPJnsIkVEpPOinOEXA/Xuvs/dTwFrgLIE7W4CHgMOJrE+ERFJkiiBPwR4PW69IbbtDDMbAswAViSvNBERSaYogW8JtnmL9fuAH7j76U89kNkCM6s2s+rTJ45GLFFERJKhR4Q2DcDQuPU84ECLNmOANWYGMBCYZmZN7r4hvpG7VwAVAL0GD2/5n4aIiKRQlMCvAoabWT7wBjALmBPfwN3zP1o2s1XAxpZhLyIi6dVm4Lt7k5ktovnqm2xgpbvvNrOFsf2atxcRyQBRzvBx903AphbbEga9u3+z82WJiEiy6U5bEZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFARAp8Mys1szozqzez8gT7y8ys1sxqzKzazK5MfqkiItIZPdpqYGbZwIPAFKABqDKzSnffE9fsj0Clu7uZFQL/CYxIRcEiItIxUc7wi4F6d9/n7qeANUBZfAN3P+7uHlvtAzgiItKlRAn8IcDrcesNsW0fY2YzzOwl4AngW4kOZGYLYlM+1adPHO1IvSIi0kFRAt8SbPvEGby7P+7uI4B/Ae5KdCB3r3D3Me4+Jrt3/3YVKiIinRMl8BuAoXHrecCB1hq7+1bgYjMb2MnaREQkiaIEfhUw3MzyzawnMAuojG9gZv/LzCy2PBroCRxKdrEiItJxbV6l4+5NZrYIeBLIBla6+24zWxjbvwL4KnCdmX0ANAIz477EFRGRLqDNwAdw903AphbbVsQtLwWWJrc0ERFJJt1pKyISCAW+iEggFPgiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBCJS4JtZqZnVmVm9mZUn2D/XzGpjP8+Z2ReSX6qIiHRGm4FvZtnAg8DVQAEw28wKWjR7BbjK3QuBu4CKZBcqIiKdE+UMvxiod/d97n4KWAOUxTdw9+fc/R+x1e1AXnLLFBGRzuoRoc0Q4PW49Qbg8k9p/2/A7ztTlIikRkVFBatXr053GZImUQLfEmzzhA3NvkRz4F/Zyv4FwAKA7H6fjViiiCTL6tWr2bq9mp6D8tNdiqRBlMBvAIbGrecBB1o2MrNC4D+Aq939UKIDuXsFsfn9XoOHJ/xPQ0RSq+egfC6YsyTdZUgH7V86vcOvjTKHXwUMN7N8M+sJzAIq4xuY2T8B64F/dfeXO1yNiIikTJtn+O7eZGaLgCeBbGClu+82s4Wx/SuAxcAA4CEzA2hy9zGpK1tERNorypQO7r4J2NRi24q45RuAG5JbmoiIJJPutBURCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFAKPBFRAKhwBcRCUSkwDezUjOrM7N6MytPsH+EmT1vZu+b2feTX6aIiHRWj7YamFk28CAwBWgAqsys0t33xDU7DHwX+JdUFCkiIp3XZuADxUC9u+8DMLM1QBlwJvDd/SBw0Mz+OSVVSpdTU1NDSUlJusuQdqqpqYF+Q9NdhqRJlMAfArwet94AXN6RNzOzBcACgOx+n+3IIaQL6FNQwnt7trB936F0lyLt1W8ofQpK0l2FpEmUwLcE27wjb+buFUAFQK/Bwzt0DEm/vkWl9C0qTXcZItJOUb60bQDifwfMAw6kphwREUmVKIFfBQw3s3wz6wnMAipTW5aIiCRbm1M67t5kZouAJ4FsYKW77zazhbH9K8zsAqAa6Ad8aGa3AAXu/m7qShcRkfaIMoePu28CNrXYtiJu+S2ap3pERKSL0p22IiKBUOCLiARCgS8iEggFvohIIBT4IiKBUOCLiARCgS8iEggFvohIIBT4IiKBUOCLiARCgS8iEggFvohIIBT4IiKBUOCLiARCgS8iEggFvohIIBT4IiKBUOCLiARCgS8iEggFvohIIBT4IiKBiBT4ZlZqZnVmVm9m5Qn2m5ktj+2vNbPRyS9VREQ6o83AN7Ns4EHgaqAAmG1mBS2aXQ0Mj/0sAP49yXWKiEgnRTnDLwbq3X2fu58C1gBlLdqUAf/Pm20HzjOzwUmuVUREOqFHhDZDgNfj1huAyyO0GQK8Gd/IzBbQ/BsAwPv7l07f1a5qM8tA4J10F5FC6l/m6s59g+7fv8939IVRAt8SbPMOtMHdK4AKADOrdvcxEd4/I6l/ma0796879w3C6F9HXxtlSqcBGBq3ngcc6EAbERFJoyiBXwUMN7N8M+sJzAIqW7SpBK6LXa0zDjjq7m+2PJCIiKRPm1M67t5kZouAJ4FsYKW77zazhbH9K4BNwDSgHjgBzIvw3hUdrjozqH+ZrTv3rzv3DdS/Vpn7J6baRUSkG9KdtiIigVDgi4gEIuWBH+GxDCVmdtTMamI/i1NdU7KY2UozO2hmCe8nyPRHTkToXyaP3VAze8bM9prZbjO7OUGbjB2/iP3L5PHLNbMdZvaXWP9+kqBNJo9flP61f/zcPWU/NH/J+zfgIqAn8BegoEWbEmBjKutIYf8mAqOBXa3snwb8nub7FMYB/5PumpPcv0weu8HA6NhyX+DlBP82M3b8IvYvk8fPgHNjyznA/wDjutH4Relfu8cv1Wf4UR7LkLHcfStw+FOaZPQjJyL0L2O5+5vu/kJs+Riwl+a7w+Nl7PhF7F/Gio3J8dhqTuyn5RUomTx+UfrXbqkO/NYeudDS+NivLr83s8tSXNPZFLX/mSzjx87MhgH/m+azqHjdYvw+pX+QweNnZtlmVgMcBP7g7t1q/CL0D9o5fqkO/CiPXHgBuNDdvwDcD2xIcU1nU6RHTmSwjB87MzsXeAy4xd3fbbk7wUsyavza6F9Gj5+7n3b3Iprv7C82s5EtmmT0+EXoX7vHL9WB3+YjF9z93Y9+dXH3TUCOmQ1McV1nS7d+5ESmj52Z5dAcho+6+/oETTJ6/NrqX6aP30fc/QiwBShtsSujx+8jrfWvI+OX6sBv87EMZnaBmVlsuThW06EU13W2dOtHTmTy2MXq/hWw191/1kqzjB2/KP3L8PH7rJmdF1s+B5gMvNSiWSaPX5v968j4RXlaZod5tMcyfA240cyagEZglse+gu7qzOy3NH9TPtDMGoAf0fzlykd968gjJ7qMCP3L2LEDrgD+FXgxNk8KcDvwT9Atxi9K/zJ5/AYDj1jzH2jKAv7T3Tda5x/50lVE6V+7x0+PVhARCYTutBURCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFA/H+1YrjP9MbgnwAAAABJRU5ErkJggg==\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/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAclElEQVR4nO3df5SUdd3/8ecbhDY3MFOpDqAsfTVa2W2yZQPsyCRom6ncnvKAY2X21ZXKzHPSr+RtnG/GQTx6yvxxf2lOmd7lJvfxVxwl9TZdTJP4YdMKaxQBxmoENyYuyIrI+/vHzA7XLrO7s7uze81c83qcw3GvuT4z8+Zy9j1vPtfnh7k7IiJS+kaEHYCIiBSGErqISEQooYuIRIQSuohIRCihi4hExFFhvfHxxx/vkyZNCuvtRURK0vr16//H3U/IdS60hD5p0iTWrVsX1tuLiJQkM3ulp3PqchERiQgldBGRiFBCFxGJiND60EUkGt555x3a2tro6OgIO5RIqaioYMKECYwaNSrv5yihi8igtLW1MWbMGCZNmoSZhR1OJLg7u3fvpq2tjaqqqryfpy4XERmUjo4OjjvuOCXzAjIzjjvuuH7/q0cJXUQGTcm88AZyTZXQRUQiQgldRGQAUqkUK1eu7Pfz4vH4kE2q1E1RibRJCx/rcrxt6edDikSiJpVKsW7dOs4555ywQ8lShS4ikfDLX/6S+vp6YrEYV1xxBX/4wx+ora2lo6ODffv2ceqpp7Jhwwaam5s544wzuOCCC6iurmbBggUcOnQIgCeffJIZM2Zw2mmnceGFF7J3714A1q5dy8yZM/n4xz9OfX09e/bsYdGiRSxfvpxYLMby5cvZt28fX/va15g2bRqf+MQn+PWvfw3A/v37mT9/PrW1tcybN4/9+/cP2TVQhS4iBXP11VeTSqUK+pqxWIzbbrut1zYvv/wyy5cv5/nnn2fUqFF84xvfYNOmTZx//vnccMMN7N+/ny996UtMnTqV5uZm1qxZQ2trKyeddBINDQ089NBDxONxFi9ezFNPPUVlZSU333wzP/zhD1m4cCHz5s1j+fLlTJs2jTfffJOjjz6aG2+8kXXr1nHnnXcCcP3113PmmWdy991388Ybb1BfX8+cOXP4yU9+wtFHH01LSwstLS2cdtppBb0+QUroIlLyfvvb37J+/XqmTZsGpKvicePGsWjRIqZNm0ZFRQW33357tn19fT2TJ08G4KKLLuK5556joqKC1tZWTj/9dAAOHDjAjBkz2LRpEx/+8Iezrz127NicMTz55JOsWLGCW2+9FUgP5/z73//Os88+y1VXXQVAbW0ttbW1Q3MRUEIXkQLqq5IeKu7OJZdcwk033dTl8R07drB3717eeecdOjo6qKysBI4cEmhmuDtnnXUWv/rVr7qca2lpyWsIobvz4IMP8tGPfvSIc8M1rDOvPnQzazCzTWa22cwW5jgfN7M9ZpbK/FlU+FBFRHKbPXs2DzzwADt37gTg9ddf55VXXqGxsZEf/OAHXHzxxVx33XXZ9mvWrGHr1q0cOnSI5cuX8+lPf5rp06fz/PPPs3nzZgDeeust/vKXvzBlyhRee+011q5dC0B7ezsHDx5kzJgxtLe3Z1/zs5/9LHfccQfuDsAf//hHAM444wzuu+8+ADZs2EBLS8uQXYc+K3QzGwncBZwFtAFrzWyFu7d2a/o7dz93CGIUEelVdXU1ixcv5uyzz+bQoUOMGjWKuXPnctRRR5FIJHj33XeZOXMmTz/9NCNGjGDGjBksXLiQl156KXuDdMSIEdxzzz1cdNFFvP322wAsXryYU045heXLl/Otb32L/fv38973vpennnqKz3zmMyxdupRYLMZ3v/tdvve973H11VdTW1uLuzNp0iQeffRRvv71r3PppZdSW1tLLBajvr5+yK6DdX6b9NjAbAbwf939s5nj7wK4+02BNnHgmv4k9Lq6OtcGFzLUNGxx6L388st87GMfCzuMvDU3N3Prrbfy6KOPhh1Kn3JdWzNb7+51udrn0+UyHtgeOG7LPNbdDDP7k5n9xsxOzfVCZtZoZuvMbN2uXbvyeGsREclXPgk9V29+97L+ReAkd/84cAfwSK4Xcveku9e5e90JJ+TcEk9EZEjF4/GSqM4HIp+E3gZMDBxPAF4LNnD3N919b+bnlcAoMzu+YFGKiEif8knoa4GTzazKzEYD84EVwQZm9iHLjMsxs/rM6+4udLAiItKzPke5uPtBM7sSeAIYCdzt7hvNbEHm/DLgi8DXzewgsB+Y733dbRURkYLKa2JRphtlZbfHlgV+vhO4s7ChiYhIf2imqIgUVPehooOVz1DTbdu2ce6557Jhw4Yujy9atIgzzjiDOXPm5HzeI488wimnnEJ1dXVBYg2bVlsUkci68cYbe0zmkE7ora3d50gOzMGDBwvyOoOhhC4ikfDuu+9y+eWXc+qpp3L22Wezf/9+vvrVr/LAAw8AsHDhQqqrq6mtreWaa67h97//PStWrODaa68lFovxt7/9jVQqxfTp06mtreWCCy7gX//6F5BePre2tpYZM2Zw7bXXMnXqVADuueceLrzwQs477zzOPvts9u7dy+zZsznttNOoqanJLqG7bds2pkyZwmWXXcbUqVO5+OKLeeqppzj99NM5+eSTWbNmTUGugRK6iETCX//6V775zW+yceNG3v/+9/Pggw9mz73++us8/PDDbNy4kZaWFm644QZmzpzJ+eefzy233EIqleIjH/kIX/nKV7j55ptpaWmhpqaG73//+wBceumlLFu2jBdeeIGRI0d2ed8XXniBe++9l6effpqKigoefvhhXnzxRZ555hm+853vZNd22bx5M9/+9rdpaWnhz3/+M01NTTz33HPceuutLFmypCDXQAldRCKhqqqKWCwGwCc/+Um2bduWPTd27FgqKiq47LLLeOihhzj66KOPeP6ePXt44403mDVrFgCXXHIJzz77LG+88Qbt7e3MnDkTgEQi0eV5Z511Fh/4wAeA9IqL119/PbW1tcyZM4dXX32Vf/7zn9n4ampqGDFiBKeeeiqzZ8/GzKipqekS62AooYtIJLznPe/J/jxy5MgufdpHHXUUa9as4Qtf+AKPPPIIDQ0Neb9uXyOwO5fkBbjvvvvYtWsX69evJ5VK8cEPfpCOjo4j4hsxYkT2eMSIEQXrf1dCF5HI27t3L3v27OGcc87htttuy+6qFFwC95hjjuHYY4/ld7/7HQC/+MUvmDVrFsceeyxjxoxh9erVANx///09vs+ePXsYN24co0aN4plnnuGVV14Z2r9YNxq2KCIFVYwrWra3tzN37lw6Ojpwd370ox8BMH/+fC6//HJuv/12HnjgAe69914WLFjAW2+9xeTJk/n5z38OwM9+9jMuv/xyKisricfjHHPMMTnf5+KLL+a8886jrq6OWCzGlClThu3vCHksnztUtHyuDActnzv0Sm353IHYu3cv73vf+wBYunQp//jHP/jxj3885O/b3+VzVaGLiPThscce46abbuLgwYOcdNJJ3HPPPWGHlJMSuohIH+bNm8e8efPCDqNPuikqIoOmtfgKbyDXVAldRAaloqKC3bt3K6kXkLuze/duKioq+vU8dbmIyKBMmDCBtrY2tK1kYVVUVDBhwoR+PUcJXUQGZdSoUVRVVYUdhqAuFxGRyFBCFxGJCCV0EZGIUB+6lJXgzFHNGpWoUYUuIhIRSugiIhGhhC4iEhHqQxcZBPXJSzFRhS4iEhFK6CIiEaGELiISEUroIiIRoYQuIhIRSugiIhGhhC4iEhF5JXQzazCzTWa22cwW9tJumpm9a2ZfLFyIIiKSjz4TupmNBO4CPgdUAxeZWXUP7W4Gnih0kCIi0rd8KvR6YLO7b3H3A8D9wNwc7b4FPAjsLGB8IiKSp3ym/o8HtgeO24BPBRuY2XjgAuBMYFrBohMZBslkkqamph7PJxIJGhsbhzEikYHJJ6Fbjse6b+99G3Cdu79rlqt55oXMGoFGgBNPPDHPEEWOFFxDBQa3jkpTUxOpVIpYLHbEuVQqBaCELiUhn4TeBkwMHE8AXuvWpg64P5PMjwfOMbOD7v5IsJG7J4EkQF1dXfcvBZFh0556nHj8FoBsMm9ubj6iXTweJ5VKEY/HAVXrUtzySehrgZPNrAp4FZgPJIIN3D275beZ3QM82j2ZixSTfa3NpN7cTiwWIxaLkUgkcrYLPq5qXYpdnwnd3Q+a2ZWkR6+MBO52941mtiBzftkQxygyJIJV+aSFj7Ekx1K4jY2N2QQej8dZvWX3Ed09IsUir/XQ3X0lsLLbYzkTubt/dfBhiRRee+px9rU2A3Bg51aYfFy/X+PAzq3saDo8FaOyOs6YWEOhQhQZFM0UlbKxr7U5nciB0eOqeuxm6UkikWD0uGzvIgd2bs1+QYgUA+1YJGVl9LgqPpRYCkBjY/9GxjQ2NrJky/jscbBSFykGSuhStnrrC9fWclKKlNAlspLJJDua7sgeH9i5tUuXSSEE+9STk1/VCBgJlfrQJbKampqyfeaQ7m6prI4X7PUrq+PZL4gDO7f2OttUZDioQpdIC/aZF9qYWEN2hIv606UYqEIXKZDOMeoapy5hUUIXEYkIJXQRkYhQH7pESnAp3FQqBWMn9v6EAtKIFwmbKnSJlM6lcCG9VkshR7X0RiNepBioQpfI6b7o1nDQiBcpBqrQRUQiQgldRCQi1OUikdDZtbJjy+4ux8VCa8PIcFCFLiISEarQpeR137ii0AtwDURwH1KA9ooabYQhQ04VupS87htXDNdQxZ5UVseJxWLZ41QqpY0wZFioQpdIGMpFuPrbHz8m1kDz0sPL9nbuRSoy1FShi4hEhBK6iEhEKKGLiESE+tClZBTb2HKRYqMKXUQkIlShS8kJjjuH4hl73hstrSvDQRW6lJzguHMojrHnvUkkElpaV4aFKnQpSUM57rzQGhsbWbJlPKCldWVoqUIXEYkIJXQRkYhQQhcRiQgldBGRiMgroZtZg5ltMrPNZnbEXR0zm2tmLWaWMrN1ZvbpwocqIiK96XOUi5mNBO4CzgLagLVmtsLdWwPNfguscHc3s1rgv4ApQxGwiIjkls+wxXpgs7tvATCz+4G5QDahu/veQPtKwAsZpEgymWRHU3pJ2lKYSNSb4OYXiURCk4ykYPLpchkPbA8ct2Ue68LMLjCzPwOPAV/L9UJm1pjpklm3a9eugcQrZaqpqamoNrEYqODmF6lUSpOMpKDyqdAtx2NHVODu/jDwsJmdAfwAmJOjTRJIAtTV1amKl34ppclEPQlufhHcok6kEPKp0NuAiYHjCcBrPTV292eBj5jZ8YOMTURE+iGfhL4WONnMqsxsNDAfWBFsYGb/y8ws8/NpwGhAe26JiAyjPrtc3P2gmV0JPAGMBO52941mtiBzfhnwBeArZvYOsB+Y5+7qUhERGUZ5Lc7l7iuBld0eWxb4+Wbg5sKGJiIi/aGZoiIiEaGELiISEVoPXYpSMpnsMkY7lUrB2Ik9P6FEBScZgSYayeCoQpei1NTUlE7iGbFYrGQnE/UkkUhkJxmBJhrJ4KlCl6IVi8Vobm7OHk9a+Fh4wQyBxsbGLtW4JhrJYKlCFxGJCFXoIkMgav+akNKgCl1EJCKU0EVEIkIJXUQkItSHLjLMgv3r25Z+PsRIJGpUoYuIRIQSuhSNZDJJPB4nHo93mVRUTjpnjsbjcZLJZNjhSIlRQpeiEZwdGovFSCQS4QY0zIIzRzVrVAZCfehSVDrGTmTb9GsBWLIFlpTReO7gzNG+Zo2qH15yUYUuIhIRSugiIhGhhC4iEhFK6CIiEaGELiISEUroIiIRoYQuodJkop5pkpH0l8ahy5DJZ6x052SiWCxGLBbjpYqa4QqvqAUnVXV+0WmvUemLErqELrjVnDaGSOvPJCORTupyERGJCCV0EZGIUJeLhG71lt3qahEpAFXoIiIRoYQuIhIRSugiIhGhhC4iEhF53RQ1swbgx8BI4KfuvrTb+YuB6zKHe4Gvu/ufChmoREMymeyyE08qlYKxE8MLqER0zhrt1F5Rw5hYQ3gBSVHqs0I3s5HAXcDngGrgIjOr7tZsKzDL3WuBHwCapyw5BbeZg/SkosrqeGjxlILg1nSQTu77WptDi0eKVz4Vej2w2d23AJjZ/cBcoLWzgbv/PtB+NTChkEFKtARnhoJmh/YlOGsU0jNHV2/ZHWJEUqzy6UMfD2wPHLdlHuvJ/wZ+k+uEmTWa2TozW7dr1678oxQRkT7lk9Atx2Oes6HZZ0gn9OtynXf3pLvXuXvdCSeckH+UIiLSp3y6XNqA4F2rCcBr3RuZWS3wU+Bz7q5/D4qIDLN8KvS1wMlmVmVmo4H5wIpgAzM7EXgI+LK7/6XwYYqISF/6rNDd/aCZXQk8QXrY4t3uvtHMFmTOLwMWAccB/2FmAAfdvW7owhaJvnzWkxcJymscuruvBFZ2e2xZ4OfLgMsKG5qIiPSHVluUIdeeepx4/BaA7O5EMjgHdm5lR9NCAJKTX9VuRgJo6r8Mg32tzdnJRLFYrMv2atJ/iUSC0eOqgHRiD868lfKmCl0KprcJQt23mVuiyUQD1tjYyJIt6akgnVW6CKhCFxGJDCV0EZGIUEIXEYkI9aHLsNC+oSJDTxW6iEhEKKGLiESEulxkSLSnHs9uwnBg59bsuGkpvOBuRolEQpOMypgSugyJfa3N2UQ+elyVdiXqQff7Cv1ds6WyOk5Nx0sA2clbSujlSwldhszocVV8KLG074YyYGNiDTQvvQOgy56jUp7Uhy4iEhFK6CIiEaGELiISEUroIiIRoYQuIhIRGuUiBZFMJtnRdEf2WGPPwxEckw4al15ulNCl33KtybKj6Y4uSVxjzwurt3VwOs+1V9QQ3AxK49LLjxK6FIzGnYcrOCYdNC69HKkPXUQkIpTQRUQiQgldRCQilNBFRCJCCV1EJCI0ykUGTGueFz+tlV5elNBlwLTmeeEVct/VRCKR/TnXmPTge/V3HXYpTkroMigae168GhsbswlcY9LLg/rQRUQiQgldRCQilNBFRCIir4RuZg1mtsnMNpvZwhznp5jZC2b2tpldU/gwRUSkL33eFDWzkcBdwFlAG7DWzFa4e2ug2evAVcC/DUWQUjyCy+RqqGJp6b60bntFDWNiDeEFJAWXT4VeD2x29y3ufgC4H5gbbODuO919LfDOEMQoRaSpqYkDO7cCWiK3lCQSCWKBtXVTqVR2DoFERz7DFscD2wPHbcCnBvJmZtYINAKceOKJA3kJKQIaqlg6Do81Hw/Tr80+3rFlYQ/tNCa9lOVToVuOx3wgb+buSXevc/e6E044YSAvISIiPcgnobcBEwPHE4DXhiYcEREZqHwS+lrgZDOrMrPRwHxgxdCGJSIi/dVnH7q7HzSzK4EngJHA3e6+0cwWZM4vM7MPAeuAscAhM7saqHb3N4cudBmI/vaVJpNJrlqszZ+j6MDOrexoSvelV1bHNeIlAvJay8XdVwIruz22LPDzDtJdMRIxnaNatPlztAT/H3aOWlJCL31anEv6pFEt0TMm1pBN4J1VupQ+Tf0XEYkIVegi0iuNUS8dqtBFRCJCFbocIZlM0tTUBGR2uhk7sfcnSMkLjnhJTn5VW9WVKFXocoSmpqbslmWxWEyjWiKusjqeHcV0YOfW7Je5lB5V6JJTLBajubkZKOw+lzK0BvL/aqAjXrq/l/rXw6cKXUQkIpTQRUQiQl0uAhx5IzS4draUF22EUbpUoQtw5I3QRCIRbkASisrqOB1jJ7J6y25Wb9nNs6vXaSOMEqIKXbKCN0KlPAVvkIKWBSg1qtBFRCJCFXqZak89Tjx+S/ZY/ebSE006Kh1K6GVqX2szqTe3Z5N4LBbjpYoajTmXLrovs9vU1KSEXsSU0MtY9z5zJXPpTsvslhb1oYuIRIQq9DKSTCbZ0ZTeTu7Azq0w+biQI5JSExyjrvHpxUcJvYwEt5MbPa5KfebSL5XVcfa1NrN6y+7M52i3EnqRUUIvM9pOTgZK/enFTwk94rS2uQyV4HBG0JDGYqCbohGntc1lKATXUAeto14sVKFHRG/7Pmptcyk0LRFQnJTQI0YzQCUswREwiURC3S8hUJdLxOxrbc52sYBWTpThUVkdzxYOqVRK3S8hUYUeAd3Hl0+fXqdVE2VYjYk10Lw0/RmMx+Oq1kOihF5Ceur/3tF0R5fx5arIJUyJRCK7nvqBnVsBuiT0nu73aI/SwVNCL1HtqcezGw90JvPO8eWNjfpFkPA0NjayZMt4IH2zVDsgDR8l9BIR7FYBeHv7BgDeM3Eqo8dVaTiiFKXK6jg1HS9lj1etWgWsyhYjGrteWEroRSw4KSj9i5BO4J3/rayOq9KRohbsW4f0Z/qqxenjt7dv4Iorrsh+xlW5D54SeoH1Nh68p3ZB7anHsxVNMIkPJoFr7LmEqevnb3y2a7Cz23D1lt2Zf3Gu6rJ/6XGpx7Of94H8LpVjH3xeCd3MGoAfAyOBn7r70m7nLXP+HOAt4Kvu/mKBY42M7t0nQW9v38AqYNasWcyaNYuXVLVIRAUnJwXvCUH69+Dt7Ruyj8VXH55boVEzPeszoZvZSOAu4CygDVhrZivcvTXQ7HPAyZk/nwL+X+a/kRbsEum0Y8vu7M/BD2FQ9+6ToPdMnMrtN3wr+4FVdS3loPvM0+4JvtOqVatYtWpVl9+74O9cUPffv3L4IsinQq8HNrv7FgAzux+YCwQT+lzgP93dgdVm9n4z+7C7/6OnF920aVOXO9+lqDMxz5o1q1/P66vy1igVKXfdE3xzpvskVxGVj1xfBFFk6RzcSwOzLwIN7n5Z5vjLwKfc/cpAm0eBpe7+XOb4t8B17r6u22s1Ap1fkR8FNhXqLzIIxwP/E3YQRULX4jBdi8N0LQ4rhmtxkrufkOtEPhW65Xis+7dAPm1w9ySQzOM9h42ZrXP3urDjKAa6FofpWhyma3FYsV+LfNZyaQOCi2hPAF4bQBsRERlC+ST0tcDJZlZlZqOB+cCKbm1WAF+xtOnAnt76z0VEpPD67HJx94NmdiXwBOlhi3e7+0YzW5A5vwxYSXrI4mbSwxYvHbqQC66ouoBCpmtxmK7FYboWhxX1tejzpqiIiJQGrYcuIhIRSugiIhGhhJ5hZteYmZvZ8WHHEhYzu8XM/mxmLWb2sJm9P+yYhpuZNZjZJjPbbGZlu1GmmU00s2fM7GUz22hm3w47prCZ2Ugz+2Nm3k1RUkIn/eElvbTB38OOJWT/DUx191rgL8B3Q45nWAWWufgcUA1cZGbV4UYVmoPAd9z9Y8B04JtlfC06fRt4OewgeqOEnvYj4P+QYzJUOXH3J939YOZwNen5BOUku8yFux8AOpe5KDvu/o/OBfbcvZ10IhsfblThMbMJwOeBn4YdS2/KPqGb2fnAq+7+p7BjKTJfA34TdhDDbDywPXDcRhknsU5mNgn4BPCHkEMJ022ki75DIcfRq7JYD93MngI+lOPUvwPXA2cPb0Th6e1auPuvM23+nfQ/ue8bztiKQF5LWJQTM3sf8CBwtbu/GXY8YTCzc4Gd7r7ezOIhh9Orskjo7j4n1+NmVgNUAX9KL+nOBOBFM6t39x3DGOKw6eladDKzS4BzgdlefpMUtIRFgJmNIp3M73P3h8KOJ0SnA+eb2TlABTDWzH7p7l8KOa4jaGJRgJltA+rcPezV1EKR2cjkh8Asd98VdjzDzcyOIn0zeDbwKullLxLuvjHUwEKQ2bTmXuB1d7865HCKRqZCv8bdzw05lJzKvg9durgTGAP8t5mlzGxZ2AENp8wN4c5lLl4G/qsck3nG6cCXgTMzn4VUpkKVIqYKXUQkIlShi4hEhBK6iEhEKKGLiESEErqISEQooYuIRIQSuohIRCihi4hExP8H2f0OWixbPc4AAAAASUVORK5CYII=\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": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAff0lEQVR4nO3df3BU9dn38fcFQlN5glqU2gEUaLE0QNzaEPnRSqxKo1V5nI4DxrZqb420WutM9Sn1ts7UOopTpz+03jdmrNW7NZV7VCxjqVqrgWqhEOw2AoqlgCVQCsVfQYiIXM8fm11ONpvkJOxmd89+XjMZ9uw53803C1x75TrfH+buiIhIdA3KdwdERCS3FOhFRCJOgV5EJOIU6EVEIk6BXkQk4o7KdwcyOf74433s2LH57oaISNFYu3btv939hEznCjLQjx07lubm5nx3Q0SkaJjZ692dU+lGRCTiFOhFRCJOgV5EJOIKskYvItHw/vvv09raSnt7e767EhllZWWMHj2aIUOGhG6jQC8iOdPa2kp5eTljx47FzPLdnaLn7uzZs4fW1lbGjRsXup1KNyKSM+3t7YwYMUJBPkvMjBEjRvT5NyQFehHJKQX57OrP+6lALyIScQr0IiJZFI/HWbZsWZ/b1dTU5GyiqG7GSsEau+C3qcdbF34xjz0RCS8ej9Pc3Mx5552X766kKKMXkUj71a9+RXV1NbFYjKuvvpo///nPVFZW0t7ezrvvvsukSZNYt24dTU1NnHHGGVx00UVUVFQwf/58Dh06BMAzzzzD9OnTOe2007j44ovZu3cvAGvWrGHGjBmceuqpVFdX8/bbb3PLLbewePFiYrEYixcv5t133+VrX/saU6dO5dOf/jS/+c1vANi/fz/z5s2jsrKSuXPnsn///py9B8roRWRAXH/99cTj8ay+ZiwW4yc/+Um351955RUWL17Miy++yJAhQ/jGN77Bxo0bufDCC7n55pvZv38/X/7yl5k8eTJNTU2sXr2aDRs2cPLJJ1NbW8vjjz9OTU0Nt912G88++yzDhg3jzjvv5Ec/+hELFixg7ty5LF68mKlTp/LOO+9w9NFHc+utt9Lc3MzPfvYzAG666SY+//nP88ADD/DWW29RXV3N2WefzX333cfRRx9NS0sLLS0tnHbaaVl9b4IU6EUksv7whz+wdu1apk6dCiSy6JEjR3LLLbcwdepUysrKuPvuu1PXV1dXM378eAAuueQSXnjhBcrKytiwYQMzZ84E4MCBA0yfPp2NGzfysY99LPXaw4cPz9iHZ555hqVLl3LXXXcBiSGn//jHP1ixYgXXXXcdAJWVlVRWVubmTUCBXopEsF4PqtkXo54y71xxdy677DLuuOOOTs/v3LmTvXv38v7779Pe3s6wYcOArkMXzQx355xzzuHXv/51p3MtLS2hhjq6O4899hif/OQnu5wbqKGnqtGLSGSdddZZPProo+zatQuAN954g9dff536+np+8IMfcOmll/Kd73wndf3q1avZsmULhw4dYvHixXz2s59l2rRpvPjii2zatAmAffv28dprrzFx4kR27NjBmjVrAGhra+PgwYOUl5fT1taWes0vfOEL3HPPPbg7AH/5y18AOOOMM3j44YcBWLduHS0tLTl7H5TRi0hkVVRUcNtttzF79mwOHTrEkCFDmDNnDkcddRR1dXV88MEHzJgxg+eee45BgwYxffp0FixYwMsvv5y6MTto0CAefPBBLrnkEt577z0AbrvtNk455RQWL17MN7/5Tfbv38+HP/xhnn32Wc4880wWLlxILBbju9/9Lt/73ve4/vrrqaysxN0ZO3YsTz75JF//+te54oorqKysJBaLUV1dnbP3wZKfMoWkqqrKtfGIpJdrglS6KQ6vvPIKn/rUp/LdjVCampq46667ePLJJ/PdlV5lel/NbK27V2W6XqUbEZGIU+lGRITEzNSampp8dyMnlNGLiEScAr2ISMSpdCNFSevgiISnjF5EJOKU0cuAUzZeunoaMtsfYf79bN26lfPPP59169Z1ev6WW27hjDPO4Oyzz87Y7oknnuCUU06hoqIiK33NJ2X0IlKSbr311m6DPCQC/YYNG7LyvQ4ePJiV1+kvBXoRibwPPviAq666ikmTJjF79mz279/P5ZdfzqOPPgrAggULqKiooLKykhtuuIE//elPLF26lBtvvJFYLMbf//534vE406ZNo7Kykosuuog333wTSCxVXFlZyfTp07nxxhuZPHkyAA8++CAXX3wxF1xwAbNnz2bv3r2cddZZnHbaaUyZMiW1XPHWrVuZOHEiV155JZMnT+bSSy/l2WefZebMmUyYMIHVq1cf8c8fKtCbWa2ZbTSzTWa2IMP5OWbWYmZxM2s2s8+GbSsikmt/+9vfuOaaa1i/fj3HHnssjz32WOrcG2+8wZIlS1i/fj0tLS3cfPPNzJgxgwsvvJAf/vCHxONxPv7xj/PVr36VO++8k5aWFqZMmcL3v/99AK644goWLVrEypUrGTx4cKfvu3LlSh566CGee+45ysrKWLJkCS+99BLPP/883/72t1Pr32zatIlvfetbtLS08Oqrr9LY2MgLL7zAXXfdxe23337EP3+vNXozGwzcC5wDtAJrzGypuwd/p/kDsNTd3cwqgf8FJoZsKxGk1SalkIwbN45YLAbAZz7zGbZu3Zo6N3z4cMrKyrjyyiv54he/yPnnn9+l/dtvv81bb73FrFmzALjsssu4+OKLeeutt2hra2PGjBkA1NXVdVpC4ZxzzuEjH/kIkFjF8qabbmLFihUMGjSI7du3869//SvVvylTpgAwadIkzjrrLMyMKVOmdOprf4XJ6KuBTe6+2d0PAI8Ac4IXuPteP7xozjDAw7YVEcm1D33oQ6nHgwcP7lQzP+qoo1i9ejVf+tKXeOKJJ6itrQ39ur2tFZZc/hjg4YcfZvfu3axdu5Z4PM5HP/pR2tvbu/Rv0KBBqeNBgwZlpb4fZtTNKGBb4LgVOD39IjO7CLgDGAkk07dQbTva1wP1ACeddFKIbkmhycaIimyPytAIH+nN3r172bdvH+eddx7Tpk3jE5/4BECn5YaPOeYYjjvuOP74xz/yuc99jl/+8pfMmjWL4447jvLyclatWsW0adN45JFHuv0+b7/9NiNHjmTIkCE8//zzvP766wPy80G4QJ9pZfwuH2PuvgRYYmZnAD8Azg7btqN9A9AAidUrQ/RLRIpMIX7YtrW1MWfOHNrb23F3fvzjHwMwb948rrrqKu6++24effRRHnroIebPn8++ffsYP348v/jFLwD4+c9/zlVXXcWwYcOoqanhmGOOyfh9Lr30Ui644AKqqqqIxWJMnDhxwH7GMIG+FRgTOB4N7OjuYndfYWYfN7Pj+9pWRCTbxo4d22kM/Q033NDlmkwjW2bOnNlleOWqVau6XDdp0qTUpiELFy6kqiqxUvDll1/O5Zdfnrru+OOPZ+XKlRn7GOzfgw8+2G3f+ytMoF8DTDCzccB2YB5QF7zAzD4B/L3jZuxpwFBgD/BWb22ltGW7VCMy0H77299yxx13cPDgQU4++eROgbpQ9Bro3f2gmV0LPA0MBh5w9/VmNr/j/CLgS8BXzex9YD8wt+PmbMa2OfpZREQG3Ny5c5k7d26+u9GjUEsguPsyYFnac4sCj+8E7gzbVkRKh7sP2CbYpaA/uwJqZqyI5ExZWRl79uzpV3CSrtydPXv2UFZW1qd2WtRMRHJm9OjRtLa2snv37nx3JTLKysoYPXp0n9oo0ItIzgwZMoRx48bluxslT6UbEZGIU0YvkaMhmyKdKaMXEYk4ZfRSMNriT/HuhqbU8bCKGspj4ReYEpHMFOglb9ID+3vbElO9PzRmMu9tW8d729alzivoi/SfAr0MqIaGBnY23gN0DuzJP5MBPfghcGDXFgAFepF+UqCXAZG8Qbqz8R4O7NrC0JHjOgX2dOWx2tTzOxsXcGDXFnY2Ht6gTBm+SHgK9NIn2dg5aujIcZxYtzD09cMqajodK8MX6RsFesm59DLM0JF9m0ATzO6ha4bfMH479fX1vb6OtjeUUqXhlZJz725oSmXhQ0eO65Kh99WwiprUh8WBXVtobGw80i6KRJoyehkQfS3X9CS9fi8iPVNGLyISccroJevSx8f3py7fF/F4nJqamsPfv2yKbtSKBCijl6wL1uQhO3X57gyrqCEWi6WO4/F4pw8ZEVFGLzmSzZp8T8pjtTQtvCd1XFNTw6rNe3L+fUWKiTJ6EZGIU0YvRS84Pn7n5j2dxthrBq2IAr1kSXANm1zffO1J8F6AZtCKJCjQS1Y0NjamAnwub772RmPsRboKFejNrBb4KTAYuN/dF6advxT4TsfhXuDr7v7XjnNbgTbgA+Cgu1dlp+tSaAbqBqyI9E2vgd7MBgP3AucArcAaM1vq7hsCl20BZrn7m2Z2LtAAnB44f6a7/zuL/ZYCkVqVUiNdRApWmIy+Gtjk7psBzOwRYA6QCvTu/qfA9auA0dnspEh/6casSLjhlaOAbYHj1o7nuvMfwO8Cxw48Y2ZrzazbJQbNrN7Mms2seffu3SG6JdKz9MXPNJFKSlWYjN4yPOcZLzQ7k0Sg/2zg6ZnuvsPMRgK/N7NX3X1Flxd0byBR8qGqqirj64v0hW7MiiSECfStwJjA8WhgR/pFZlYJ3A+c6+6pgq277+j4c5eZLSFRCuoS6KVwpa/jnnSk68znW/Dn0tr0EmVhSjdrgAlmNs7MhgLzgKXBC8zsJOBx4Cvu/lrg+WFmVp58DMwG1mWr85Jf2V5nXkRyo9eM3t0Pmtm1wNMkhlc+4O7rzWx+x/lFwC3ACOC/zAwOD6P8KLCk47mjgEZ3fyonP4nkRTENqdS+s1KqQo2jd/dlwLK05xYFHl8JXJmh3Wbg1CPso8gR076zUso0M1ZKQqZ9Z0VKhVavFBGJOGX0ImgEjkSbAr2ENtBbBIpIdijQS2jJ4ZTJ4F7sQyq1PIKUCgV66ZNiGk7ZE61bL6VEgV5KkpZHkFKiUTciIhGnQC8iEnEK9CIiEacavQidR+A0jN9OfX23WyeIFB1l9NKjhoYGdjYuYGfjgtTolKhJ36CksbExzz0SyS5l9NKjxsbG1Nj5Yh833x2NwJGoU6CXXkVl7LxIqVLpRkQk4pTRi6RZtXmPFjmTSFGgl4ySgW7n5j29XBk96TtRaRSOFDuVbkQCgiNwQKNwJBqU0YsEaCcqiSJl9CIiEaeMXrpITJK6B9DmIiJRoEAvQOet9HY23hP5SVIipSRUoDezWuCnwGDgfndfmHb+UuA7HYd7ga+7+1/DtJXCpElSmWnYpRSjXgO9mQ0G7gXOAVqBNWa21N03BC7bAsxy9zfN7FygATg9ZFuRghaPx6mpqQGgrWyKdqKSohPmZmw1sMndN7v7AeARYE7wAnf/k7u/2XG4Chgdtq1IIRtWUUMsFgMSAT+4ObpIsQgT6EcB2wLHrR3Pdec/gN/1ta2Z1ZtZs5k17969O0S3RHKvPFZLU1MTTU1NqYAvUmzCBHrL8JxnvNDsTBKBPlmvD93W3Rvcvcrdq0444YQQ3RIRkTDC3IxtBcYEjkcDO9IvMrNK4H7gXHff05e2IiKSO2Ey+jXABDMbZ2ZDgXnA0uAFZnYS8DjwFXd/rS9tRUQkt3rN6N39oJldCzxNYojkA+6+3szmd5xfBNwCjAD+y8wADnaUYTK2zdHPIpJz2nJQilGocfTuvgxYlvbcosDjK4Erw7aVwhKcCQuaDduduro6VnWs5plc7EyBXoqBZsZKp+0CAc2G7UZ9fT23b04MGtNiZ1JMFOgF0ExYkSjT6pUiIhGnQC8iEnEK9CIiEacavUg/BRc7g8SoHI3CkUKkQC/Si+DSxEnDKmqY0v5y6jgejwMo0EtBUqAX6YfyWC1NCw/PPQhm9iKFRjV6EZGIU0ZfohoaGmhsbAQ6yg7Dx/TcQLrotP3i5j1MGz8ij70R6Z4y+hLV2NiYqivHYjHNhBWJMGX0JSwWi9HU1ARkvuEofRMchfNyYMtB7S0r+aZAX8JWbd6jAJ8lwVE48Xic9uF7tLesFAyVbkSyQFsOSiFToBcRiTgFehGRiFOgFxGJON2MFckBbTkohUQZvUiW1dXVpXbrSm45KJJPyuhLiGbDDgxtOSiFRhl9CdFsWJHSpIy+xGg2rEjpCZXRm1mtmW00s01m1uV3UTObaGYrzew9M7sh7dxWM3vZzOJm1pytjouISDi9BnozGwzcC5wLVACXmFlF2mVvANcBd3XzMme6e8zdq46ksyLFaMWqZspOmpL6amhoyHeXpMSEKd1UA5vcfTOAmT0CzAE2JC9w913ALjPT6k1SsrrbiSooOQpHwy1lIIUJ9KOAbYHjVuD0PnwPB54xMwfuc3elM1IyymO1nRY30ygcyYcwgd4yPOd9+B4z3X2HmY0Efm9mr7r7ii7fxKweqAc46aST+vDyIiLSkzA3Y1uB4IDr0cCOsN/A3Xd0/LkLWEKiFJTpugZ3r3L3qhNOOCHsy4uISC/CBPo1wAQzG2dmQ4F5wNIwL25mw8ysPPkYmA2s629nRUSk73ot3bj7QTO7FngaGAw84O7rzWx+x/lFZnYi0AwMBw6Z2fUkRugcDywxs+T3anT3p3Lyk0gXwZmwkJgNq7XS8y+4E1VdXZ1uzErOhZow5e7LgGVpzy0KPN5JoqST7h3g1CPpoPRfciZsMrjHYjHq6ury26kSl74TFaBALzmnmbERF5wJK/lXHqulaeE9AKmsXiTXFOhLjJY9yL/k38HOzXuYNn5EnnsjpUCLmomIRJwCvUgeJW/M1tTUaGkEyRmVbkTyRDdmZaAo0IvkiW7MykBR6UZEJOIU6EVEIk6lm4hJ3xdWM2ELW3CoZfAYYOtCrfot2aFAHzHB2bCaCVtcDuza0mkZ44bx23VzVrJCgT6CNBu2+GiDEsklBXqRAqANSiSXdDNWRCTiFOhFRCJOgV6kQGl5BMkW1ehFCpCWR5BsUqAXKUBaHkGySaUbEZGIU6CPgIaGhlQtN/lrvkSL6vVyJFS6iYD02bAvl03RTlIREpzdrHq99IcCfUQEZ8MqyEfD4b/HUWzt+LtVvV76Q6UbEZGICxXozazWzDaa2SYz6zI328wmmtlKM3vPzG7oS1sREcmtXks3ZjYYuBc4B2gF1pjZUnffELjsDeA64P/2o62I9EHyxmxSXV2davbSozAZfTWwyd03u/sB4BFgTvACd9/l7muA9/vaVkTCq6ur67THQDweT+0/INKdMDdjRwHbAsetwOkhX/9I2opImvr6+k7Zu27OShhhMnrL8JyHfP3Qbc2s3syazax59+7dIV9eRER6EyajbwXGBI5HAztCvn7otu7eADQAVFVVhf0gKUnB7QIh8et7+/AxGlYZYT1tMRis2ateL5mEyejXABPMbJyZDQXmAUtDvv6RtJVuJCdIJcVisS47FElpCNbsVa+X7vSa0bv7QTO7FngaGAw84O7rzWx+x/lFZnYi0AwMBw6Z2fVAhbu/k6ltjn6WkpK+XaCy+dLR+e9ak6mkd6Fmxrr7MmBZ2nOLAo93kijLhGorIiIDRzNjRSJEi59JJlrrRiQitPiZdEeBXiQigmPsVa+XIJVuREQiThm9SERpTRxJUqAvEsFJUslNRkS6E6zXg2r2pU6Bvkik7yKV/h9ZJEhr4kiQAn0RSZ8kJdIXWiqhdCnQi0RIdzOk28qmkKz2xeNxVm3ew+2bR6XOp6+fI9GiUTciJaA8VktTUxNNTU26v1OClNEXKa1tIyJhKdCLFLn+fOgf2LWFnY2Ht3BuGL9dNfsIU+lGpMTU1dUxdOS41PGBXVs6LW88dsFvU18SDcroC5jGzks2HQ7cozixbmHq+WBmL9GkjL6ABTcY0dh5yaXgqpdt8afy3R3JMmX0BU5j5yXXhlXUMKX9ZSC5LeUeymO1ee6VZJMCvUiJK4/V0rTwHiAxg3bFquZUOUc3aaNBpRsRSQneqE2/SSvFSxm9iKTU19enZszqJm10KKMXkW6tWNVM2UlTUl/anrA4KdAXkIaGhtTIh5qamtSIG5F8GFZR0+N4eykeKt0UkOBSxKAhlZJf5bHaTqNvdjYu0AqYRUqBvsBoOKUUqvRhmKCNTIpFqEBvZrXAT4HBwP3uvjDtvHWcPw/YB1zu7i91nNsKtAEfAAfdvSprvReRrAiz3EH6MEwpHr3W6M1sMHAvcC5QAVxiZhVpl50LTOj4qgf+O+38me4eU5AXiY7gbFrdpC1sYTL6amCTu28GMLNHgDnAhsA1c4D/cXcHVpnZsWb2MXf/Z9Z7LCJ5F7x3pDJO4Qsz6mYUsC1w3NrxXNhrHHjGzNaaWbf/Esys3syazax59+7dIboVDcGRNhplI8Wivr6+00YmwexeGX7hCRPoLcNz3odrZrr7aSTKO9eY2RmZvom7N7h7lbtXnXDCCSG6FQ1auEyKXV1dXaeVVePxuIZhFpgwpZtWYEzgeDSwI+w17p78c5eZLSFRClrR3w5HUdiRNlofXApRfX19p7JN8rdTDcMsHGEC/RpggpmNA7YD84D0tHMpcG1H/f504G13/6eZDQMGuXtbx+PZwK3Z676I5Esw8QhuLq76feHpNdC7+0EzuxZ4msTwygfcfb2Zze84vwhYRmJo5SYSwyuv6Gj+UWBJYvQlRwGN7q7FrkUiLJjhK7svDKHG0bv7MhLBPPjcosBjB67J0G4zcOoR9lFECkRfy4fB7H758uUsX768U/1egX9gaGbsAAtuDwjaIlCiLZjdZ/q3n7xGckuBfoBpPRspBZkz/1Ew7cbUUfvmBanVMQHuvvmbCvo5okCfB1rPRqKmPyPChlXUpB6/t20dV199dSrjV0knuxToC5yGVEpUBVfHbIs/lVowLVnLv+62e1LXKts/MlqPXkTyrjxWm5ppe9999/GhMZNT57QO/pFTRj8AgjehdPNVJLPDv72O4sS6wwvkah38I6dAPwCCN2B7u/mqUo1IZ8F18NOHaKYH/e4mcZU6BfoBohuwIv0TXAc/+NtxpnH5bWVTOu2KJQkK9CJSNHoal798+XJgOe9uaAJgRPypVNAv9exegT5H+lKXV7lGpO/SF1NraGhIjdR5b9s63tu2LhX0G8Zv77aun/7/L4ofCgr0OdKXuryI9Ky72nvnIH34Jm5b/KlUkE8fo5+uFMo9CvQ5pLq8SH50N0Y/XXq5B3rO/ouVAn2WaA0bkYHR11Jn8GZuumC5B3rO/ot5hI8CfZZoDRuR4lNfX8/tmw/vjNpd9l/sI3wU6LNIpRqRwhT2t4DyWC1bORy8k5l6byN8alb9MHWuECd0KdAfAc14FYm24GzdrYEkLr3kk5Qp8w/q6UMgl6UgBfojoJE1IqWjuxE+AE3dZP5BvX0I7Ny8J/U42zeEFej7KFMWr3KNSGnrLvMP6ulDICj9hvCqwAcAwLTxI1KPw5aJFOj7SFm8iPRH+o3fdCdOS/zZ03DQoN5+QwhSoO9Fd8MmlcWLSCZHOtM2fTho+uuFKROlU6DPIH3hJIBZs2YBfRs2GW42n4hEWa7+v6cvAWFm3V6rQN+hu+A+a9asrAyXUnAXkWwLG1dCBXozqwV+CgwG7nf3hWnnreP8ecA+4HJ3fylM23zJPC6278G9mGbHiUhx6CmA9ydp7DXQm9lg4F7gHKAVWGNmS919Q+Cyc4EJHV+nA/8NnB6y7RHrS60qKb0kk63MXUSk0ITJ6KuBTe6+GcDMHgHmAMFgPQf4H3d3YJWZHWtmHwPGhmjbxcaNG1PbhoWRHrTDOJLA3t0nqsozIlKIwgT6UcC2wHEriay9t2tGhWwLgJnVA8mou3f58uUbQ/Stk2TA78v1V199dXenjwf+3dc+RJTei8P0Xhym9+KwQngvTu7uRJhAn+lWroe8JkzbxJPuDUBDiP4MCDNrdveqfPejEOi9OEzvxWF6Lw4r9PciTKBvBcYEjkcDO0JeMzREWxERyaFBIa5ZA0wws3FmNhSYByxNu2Yp8FVLmAa87e7/DNlWRERyqNeM3t0Pmtm1wNMkhkg+4O7rzWx+x/lFwDISQys3kRheeUVPbXPyk2RfwZSRCoDei8P0Xhym9+Kwgn4vLDFQRkREoipM6UZERIqYAr2ISMQp0IdgZjeYmZvZ8fnuS76Y2Q/N7FUzazGzJWZ2bL77NJDMrNbMNprZJjNbkO/+5IuZjTGz583sFTNbb2bfynef8s3MBpvZX8zsyXz3pTsK9L0wszEklnD4R777kme/Bya7eyXwGvDdPPdnwASW8jgXqAAuMbOK/PYqbw4C33b3TwHTgGtK+L1I+hbwSr470RMF+t79GPh/dDPRq1S4+zPufrDjcBWJORGlIrUMiLsfAJJLeZQcd/9ncsFCd28jEeC6300j4sxsNPBF4P5896UnCvQ9MLMLge3u/td896XAfA34Xb47MYC6W+KjpJnZWODTwJ/z3JV8+gmJRPBQnvvRo5Jfj97MngVOzHDqP4GbgNkD26P86em9cPffdFzznyR+fX94IPuWZ6GX8igVZvZ/gMeA6939nXz3Jx/M7Hxgl7uvNbOaPHenRyUf6N397EzPm9kUYBzw146dW0YDL5lZtbvvHMAuDpju3oskM7sMOB84y0trAkaYZUBKhpkNIRHkH3b3x/PdnzyaCVxoZucBZcBwM/uVu385z/3qQhOmQjKzrUCVu+d7hbq86NhA5kfALHffne/+DCQzO4rEDeizgO0klvaoK6JZ3lnTscnQQ8Ab7n59nrtTMDoy+hvc/fw8dyUj1eglrJ8B5cDvzSxuZovy3aGB0nETOrmUxyvA/5ZikO8wE/gK8PmOfwfxjoxWCpgyehGRiFNGLyIScQr0IiIRp0AvIhJxCvQiIhGnQC8iEnEK9CIiEadALyIScf8fSoP3/56JUvMAAAAASUVORK5CYII=\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\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_{-\\infty}^r te^{-t^2/2}dt = -e^{-r^2/2}\n",
|
|
"$$\n",
|
|
"which can be calculated using a substitution with $z := t^2$.\n",
|
|
"\n",
|
|
"Now the inversion.\n",
|
|
"\n",
|
|
"$$\n",
|
|
"p := F_R(r) = -e^{-r^2/2} \\implies r = \\sqrt{\\ln{p^2}}\n",
|
|
"$$\n",
|
|
"for $p \\in [0, 1]$."
|
|
]
|
|
},
|
|
{
|
|
"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/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAglElEQVR4nO3df3CcVfn38ffV0hiJKWJpwekPm2oRQxvWkoa2KCwCNSDQh/FhWoqC+IVY5YfMCA8VkXkERsvIqIB8v3UHER5lbR2g2IEqiJCWX7FNdQ1toVjTYkOpKUGghYZQej1/7GbZpJvkTpPN7t75vGYy2fu+z9k9exeuPTl7znXM3RERkfAake8GiIhIbinQi4iEnAK9iEjIKdCLiIScAr2ISMgdku8GZHPEEUf45MmT890MEZGisX79+tfcfWy2awUZ6CdPnkxjY2O+myEiUjTM7OWermnoRkQk5BToRURCToFeRCTkCnKMXkTC4b333qOlpYX29vZ8NyU0SktLmTBhAqNGjQpcR4FeRHKmpaWF8vJyJk+ejJnluzlFz91pa2ujpaWFioqKwPU0dCMiOdPe3s6YMWMU5AeJmTFmzJh+/4WkQC8iOaUgP7gO5n4GCvRmVmtmm81si5kt7qXcTDN738z+d3/riohIbvQZ6M1sJHAncAZQCZxvZpU9lLsFeLS/dUVEwiKRSLBq1ap+14tGozlbKBqkR18DbHH3ZnfvAJYB87KUuwJ4AGg9iLoiRWHy4kfSPyLZHGygz6UggX48sD3juCV1Ls3MxgPnAkv7WzfjOerMrNHMGnft2hWgWSK5F4vFiEaj6Z+d8cXpn1gslu/mSQC/+c1vqKmpIRKJ8I1vfIO//OUvVFVV0d7ezttvv82xxx7Lhg0bqK+v56STTuLcc8+lsrKSRYsWsX//fgAee+wxZs+ezYwZMzjvvPPYs2cPAOvWrWPOnDkcd9xx1NTU8Oabb3LDDTewfPlyIpEIy5cv5+233+brX/86M2fO5LOf/Sy///3vAdi7dy8LFiygqqqK+fPns3fv3pzdgyDTK7ON/Hfff/BnwLXu/n63LwqC1E2edI8BMYDq6mrtbygFIR6Ps6ahkZJxXaeydbRuJR6PU1dXl6eWFZ+rrrqKRCIxqM8ZiUT42c9+1uP1F154geXLl/PMM88watQovvWtb7F582bOOeccrr/+evbu3ctXvvIVpk2bRn19PWvXrmXTpk184hOfoLa2lgcffJBoNMrNN9/M448/TllZGbfccgs/+clPWLx4MfPnz2f58uXMnDmTt956i0MPPZQbb7yRxsZGfv7znwNw3XXX8YUvfIG7776bN954g5qaGk477TR+8YtfcOihh9LU1ERTUxMzZswY1HuTKUigbwEmZhxPAHZ0K1MNLEsF+SOAM81sX8C6IgUjFosRj8fTx4lEgpJxFRy1cEmXcjvji1nT0EjppOnpc7dff4UCf4H585//zPr165k5cyaQ7EWPGzeOG264gZkzZ1JaWsrtt9+eLl9TU8OUKVMAOP/883n66acpLS1l06ZNnHjiiQB0dHQwe/ZsNm/ezMc//vH0c48ePTprGx577DFWrlzJrbfeCiSnnP7rX/9izZo1XHnllQBUVVVRVVWVm5tAsEC/DphqZhXAK8ACYGFmAXdPd3fM7B7gYXd/yMwO6auuSCGJx+MkEgkikQiQ7DE+Xzr9gHJlldEux+rh9623nneuuDsXXXQRP/rRj7qc37lzJ3v27OG9996jvb2dsrIy4MCpi2aGu3P66afz29/+tsu1pqamQFMd3Z0HHniAT3/60wdcG6qpp32O0bv7PuBykrNpXgB+5+4bzWyRmS06mLoDb7ZI7rSPnsi2Wdekf8ojtQeUKY/UctTCJemf7kM7UhhOPfVU7r//flpbk3NEXn/9dV5++WXq6uq46aabuOCCC7j22mvT5deuXcvWrVvZv38/y5cv53Of+xyzZs3imWeeYcuWLQC88847vPTSSxxzzDHs2LGDdevWAbB792727dtHeXk5u3fvTj/nF7/4Re644w7ckyPSf/vb3wA46aSTuO+++wDYsGEDTU1NObsPgVIguPsqYFW3c92/eO08/7W+6oqIDIXKykpuvvlm5s6dy/79+xk1ahTz5s3jkEMOYeHChbz//vvMmTOHJ554ghEjRjB79mwWL17M888/n/5idsSIEdxzzz2cf/75vPvuuwDcfPPNHH300SxfvpwrrriCvXv38uEPf5jHH3+cU045hSVLlhCJRPjud7/L97//fa666iqqqqpwdyZPnszDDz/MN7/5TS6++GKqqqqIRCLU1NTk7D5Y56dMIamurnZtPCJDJXNcPpFI0D564gFj8n3ZGV9M6Vvb00M+Cxcu1DAOyS9DP/OZz+S7GYHU19dz66238vDDD+e7KX3Kdl/NbL27V2crr6RmMux1mVkzeuIB4+9BlFVGeXtTPQ3NbXS0bgVQoJeCoUAvAlln1vRHeaQ2PZa/M65MH8Woc61EGCnQy7DVubp1Z3NbnlsiklvKXikiEnIK9CI50NDcppw4UjA0dCPDUiwWY2f8DiC52Gmw58F3tG5Nj9XHpryiL2YlrxToZViKx+PpAF8yruKgZtr0JPO5tGK2q8H+C2fbki/1XWbbNs466yw2bNjQ5fwNN9zASSedxGmnnZa13kMPPcTRRx9NZWXxZ1ZXoJdha6AzbXqiGTjF4cYbb+z1+kMPPcRZZ501KIF+3759HHJI/sKtxuhFJPTef/99Lr30Uo499ljmzp3L3r17+drXvsb9998PwOLFi6msrKSqqoqrr76aZ599lpUrV3LNNdcQiUT45z//SSKRYNasWVRVVXHuuefyn//8B0imKq6qqmL27Nlcc801TJs2DYB77rmH8847j7PPPpu5c+eyZ88eTj31VGbMmMH06dPT6Yq3bdvGMcccwyWXXMK0adO44IILePzxxznxxBOZOnUqa9euHfD7V6AXkdD7xz/+wWWXXcbGjRv56Ec/ygMPPJC+9vrrr7NixQo2btxIU1MT119/PXPmzOGcc87hxz/+MYlEgk9+8pNceOGF3HLLLTQ1NTF9+nR+8IMfAHDxxRezdOlSnnvuOUaOHNnldZ977jnuvfdennjiCUpLS1mxYgV//etfefLJJ/nOd76Tzn+zZcsWvv3tb9PU1MSLL75IPB7n6aef5tZbb+WHP/zhgN+/Ar2IhF5FRUU6PcXxxx/Ptm3b0tdGjx5NaWkpl1xyCQ8++CCHHnroAfXffPNN3njjDU4++WQALrroItasWcMbb7zB7t27mTNnDpBMfZHp9NNP52Mf+xiQzGJ53XXXUVVVxWmnncYrr7zCv//973T7pk+fzogRIzj22GM59dRTMTOmT5/epa0HS4FeRELvQx/6UPrxyJEj2bdvX/r4kEMOYe3atXz5y1/moYceorb2wGylPekrV1hn+mOA++67j127drF+/XoSiQRHHnkk7e3tB7RvxIgR6eMRI0Z0aevB0pexMixk21CE0RN7rjCIEolEl6X1SnhWWPbs2cM777zDmWeeyaxZs/jUpz4F0CXd8GGHHcbhhx/OU089xec//3l+/etfc/LJJ3P44YdTXl5OQ0MDs2bNYtmyZT2+zptvvsm4ceMYNWoUTz75JC+//PKQvD9QoJdh4oAtAQ8yeVl/lVVGmd7+fPq4cyu94Rrog0yHHGq7d+9m3rx5tLe34+789Kc/BWDBggVceuml3H777dx///3ce++9LFq0iHfeeYcpU6bwq1/9CoBf/vKXXHrppZSVlRGNRjnssMOyvs4FF1zA2WefTXV1NZFIhGOOOWbI3qPSFMuwEI1GaWhuy8l0yr5kBrfOnn19ff2QtyMfiilN8cHas2cPH/nIRwBYsmQJr776KrfddltOX7O/aYoDjdGbWa2ZbTazLWZ2wMRgM5tnZk1mljCzRjP7XMa1bWb2fOe1fr4fEZGC9sgjjxCJRJg2bRpPPfUU119/fb6bdIA+h27MbCRwJ3A6yc2+15nZSnfflFHsz8BKd3czqwJ+B2T+XXKKu782iO0WESkI8+fPZ/78+fluRq+C9OhrgC3u3uzuHcAyYF5mAXff4x+MAZUBhTceJCJ5UYjDw8XsYO5nkC9jxwPbM45bgBO6FzKzc4EfAeOAzG9cHHjMzBz4hbvHsr2ImdUBdQCTJk0K1HiRYpCZ32VncxuzpozJY2uGVmlpKW1tbYwZMwYzy3dzip6709bWRmlpab/qBQn02f51DvhIcfcVwAozOwm4CejMFHSiu+8ws3HAn8zsRXdfk6V+DIhB8svYoG9ARArXhAkTaGlpYdeuXfluSmiUlpYyYcKEftUJEuhbgMwJxxOAHT0Vdvc1ZvZJMzvC3V9z9x2p861mtoLkUNABgV5EwmfUqFFUVAxuCmjpvyBj9OuAqWZWYWYlwAJgZWYBM/uUpf4uM7MZQAnQZmZlZlaeOl8GzAW65goVyZFYLJbeB7Rz/nohWNPQSOmk6ZROmk4slnUkU2RQ9Rno3X0fcDnwKPAC8Dt332hmi8xsUarYl4ENZpYgOUNnfurL2SOBp83s78Ba4BF3/2MO3ofIAeLxeDrARyKRIVkg1Zeyymh60VZnrnqRXNOCKQmt7ouTCm1bv53xxcyaMmbYLJ6S3OptwZRSIEiode7dKjKcKXuliEjIKdCLiIScAr2ISMgp0IuIhJwCvYhIyCnQi4iEnAK9hEqhrobtSec2g9FoVKtkJWcU6CVUCnE1bE/KKqO0j55IQ3MbaxoatUpWckYLpiR0IpFIwa6GzVQeqaU8UgskV8mK5Ip69CIiIadALyIScgr0IiIhp0AvIhJy+jJWQqdYM1Z2b/e2JV/qpbRIcOrRi4iEnAK9iEjIBRq6MbNa4DZgJHCXuy/pdn0ecBOwH9gHXOXuTwepKzJQmcMdO5vb8tgSkcLUZ4/ezEaS3Af2DKASON/MKrsV+zNwnLtHgK8Dd/WjroiI5FCQHn0NsMXdmwHMbBkwD9jUWcDd92SULwM8aF2RgYjFYuyM35E+7mjdmt58u9h0tG7tskI2NuUV6urq8tgiCYsgY/Tjge0Zxy2pc12Y2blm9iLwCMlefeC6qfp1ZtZoZo27du0K0nYR4vE4Ha1b08cl4yoKOr9NT8oqo10+oDpatyr3jQyaID16y3LODzjhvgJYYWYnkRyvPy1o3VT9GBADqK6uzlpGJJuScRUctbC4v/rJzHsDyn0jgytIj74FmJhxPAHY0VNhd18DfNLMjuhvXRERGXxBAv06YKqZVZhZCbAAWJlZwMw+ZWaWejwDKAHagtQVEZHc6nPoxt33mdnlwKMkp0je7e4bzWxR6vpS4MvAhWb2HrAXmO/uDmStm6P3IiIiWQSaR+/uq4BV3c4tzXh8C3BL0LoiIjJ0tDJWRCTkFOhFClRnkrNiTNAmhUWBXkQk5BToRURCToFeRCTktPGIFJ1YLJZOD5BIJGD0xN4rFKnM3DfKeyMDoR69FJ14PJ4M8EAkEinK3DZ9ycx9o7w3MlDq0UtRikQi1NfXA4RyVkpm7hvlvZGBUo9eRCTkFOhFREJOQzdSlDoXE4lI39SjFxEJOQV6EZGQ09CNFIXMYZqdzW15bIlI8VGPXqQIKMGZDIQCvYhIyCnQi4iEXKAxejOrBW4juR3gXe6+pNv1C4BrU4d7gG+6+99T17YBu4H3gX3uXj04TZfhIhaLsTN+R/q4o3VrOj3AcKG8NzIQffbozWwkcCdwBlAJnG9mld2KbQVOdvcq4CYg1u36Ke4eUZCXgxGPx+lo3Zo+LhlXEcr8Nj1R3hsZqCA9+hpgi7s3A5jZMmAesKmzgLs/m1G+AZgwmI0UKRlXwVELl/RdMISU90YGKsgY/Xhge8ZxS+pcT/4L+EPGsQOPmdl6M+vx700zqzOzRjNr3LVrV4BmiYhIEEF69JblnGctaHYKyUD/uYzTJ7r7DjMbB/zJzF509zUHPKF7jNSQT3V1ddbnFxGR/gvSo28BMnd2mADs6F7IzKqAu4B57p5e0eLuO1K/W4EVJIeCRERkiAQJ9OuAqWZWYWYlwAJgZWYBM5sEPAh81d1fyjhfZmblnY+BucCGwWq8iIj0rc+hG3ffZ2aXA4+SnF55t7tvNLNFqetLgRuAMcB/mxl8MI3ySGBF6twhQNzd/5iTdyIiIlkFmkfv7quAVd3OLc14fAlwSZZ6zcBxA2yjiIgMgJKaScHqzOuiJGYiA6MUCCIiIadALyISchq6kYKUmd9mOOa26U0ikSAajaaPFy5cqNw30iv16KUgZea3GW65bXpTVhklEomkjxOJhHLfSJ/Uo5eCNZzz2/SkPFJL/ZIPMnlm9uxFeqIevYhIyKlHL1Jkuu+fO2vKmDy2RoqBevQiIiGnQC8iEnIK9CIiIadALyIScgr0IiIhp0AvIhJyCvQiIiGnQC8iEnIK9FIwYrEY0WiUaDRKIpHId3OKRmeSs2g0SiwWy3dzpAAFCvRmVmtmm81si5ktznL9AjNrSv08a2bHBa0r0ikej6cDfCQSUSKzAMoqo7SPnkhDcxtrGhqV4Eyy6jMFgpmNBO4ETgdagHVmttLdN2UU2wqc7O7/MbMzgBhwQsC6ImmRSIT6+nqg61J/ya48Ukt5pBaAnXH1oyS7ILluaoAtqf1fMbNlwDwgHazd/dmM8g3AhKB1RTI1NLcpwIsMsiBDN+OB7RnHLalzPfkv4A/9rWtmdWbWaGaNu3btCtAsEREJIkigtyznPGtBs1NIBvpr+1vX3WPuXu3u1WPHjg3QLBERCSLI0E0LMDHjeAKwo3shM6sC7gLOcPe2/tQVEZHcCdKjXwdMNbMKMysBFgArMwuY2STgQeCr7v5Sf+qKiEhu9dmjd/d9ZnY58CgwErjb3Tea2aLU9aXADcAY4L/NDGBfahgma90cvRcREcki0A5T7r4KWNXt3NKMx5cAlwStKyIiQ0crY0VEQk6BXiREOtchaC2CZNLm4JI3sVisy5L9RCIBoyf2XEH61NG6Nb1CNjblFerq6vLcIikE6tFL3sTjcdY0NNLQ3EZDcxvtoycqv80AlFVGKRlXASQDvvLeSCf16CWvSsZVcNTCJfluRigo7430RD16EZGQU6AXEQk5BXoRkZBToBcRCTkFehGRkFOgFxEJOQV6EZGQU6AXEQk5LZiSIZWZ9kApD3IrkUgQjUbTxwsXLlRKhGFKgV6GVGfag5JxFaCUBzlTVhllevvz6eNEIgGgQD9MKdDLkFPag9wrj9RSv+SO9HFmz16Gn0CB3sxqgdtI7hJ1l7sv6Xb9GOBXwAzge+5+a8a1bcBu4H1SO08NTtNFpDeZqYp3Nrcxa8qYPLZG8qnPQG9mI4E7gdNJbva9zsxWuvumjGKvA1cC/6uHpznF3V8bYFtFROQgBJl1UwNscfdmd+8AlgHzMgu4e6u7rwPey0EbRURkAIIE+vHA9ozjltS5oBx4zMzWm5m+CRIRGWJBxugtyznvx2uc6O47zGwc8Ccze9Hd1xzwIskPgTqASZMm9ePpRUSkN0F69C1A5mTnCcCOoC/g7jtSv1uBFSSHgrKVi7l7tbtXjx07NujTi4hIH4IE+nXAVDOrMLMSYAGwMsiTm1mZmZV3PgbmAhsOtrEiItJ/fQ7duPs+M7sceJTk9Mq73X2jmS1KXV9qZkcBjcBoYL+ZXQVUAkcAK8ys87Xi7v7HnLwTERHJKtA8endfBazqdm5pxuOdJId0unsLOG4gDZTip7QHhWFNQyOlk6YDcPv1V2iV7DCipGaSc/F4PL0EPxKJKO1BHpRVRpNpJ4CO1q3pD14ZHpQCQYZE++iJbJt1DQDleW7LcFQeqaU8UgvAzvjiPLdGhpp69CIiIadALyIScgr0IiIhp0AvIhJyCvQiIiGnQC8iEnIK9CIiIad59DLoYrEYV978wTZ2Ha1b04t1pDBo4/DhRYFeBl08Hu8S3EvGVWg1bAEpq4zy9qZ6GprbgOQHMWjj8DBToJec0AbghStzlSxopexwoDF6EZGQU6AXEQk5BXoRkZBToBcRCTkFehGRkAsU6M2s1sw2m9kWMzvgK3ozO8bMnjOzd83s6v7UFRGR3Ooz0JvZSOBO4AyS+8Ceb2aV3Yq9DlwJ3HoQdSUEYrEY0WiUaDSa3k1KikfnAqpoNEosFst3c2SQBenR1wBb3L3Z3TuAZcC8zALu3uru64D3+ltXwkHbBRavssoo7aMn0tDcxpqGRm0zGEJBFkyNB7ZnHLcAJwR8/sB1zawOqAOYNGlSwKeXQqLtAouTthkMvyA9estyzgM+f+C67h5z92p3rx47dmzApxcRkb4ECfQtwMSM4wnAjoDPP5C6IiIyCIIE+nXAVDOrMLMSYAGwMuDzD6SuiIgMgj7H6N19n5ldDjwKjATudveNZrYodX2pmR0FNAKjgf1mdhVQ6e5vZaubo/ciIiJZBMpe6e6rgFXdzi3NeLyT5LBMoLoiIjJ0tDJWDprmzoeT5tSHjwK9HDTNnQ+fssookUgESAZ8zakPB208IgOiufPhUh6ppX5JchvIzK0GpbipRy8iEnIK9CIiIaehGxHpYvLiRwDY2dzGrClj8twaGQzq0YuIhJwCvYhIyGnoRgKLxWJdptslEgkYPbHnClL01jQ0Ujppevr49uuvoK6uLo8tkoOhHr0EljlvHjR3PuzKKqOUjKtIH3e0btW8+iKlHr30S+a8edDc+TDLzFMPylVfzNSjFxEJOQV6EZGQU6AXEQk5BXrplTJUSiZltixO+jJWehWPx1nT0JicfTF6ombZDGNllVHe3lRPQ3MbHa1bATTVskgo0EufSsZVcNTCJfluhuRZ5iwczcApLoGGbsys1sw2m9kWMzvgX9iSbk9dbzKzGRnXtpnZ82aWMLPGwWy8iIj0rc8evZmNBO4ETgdagHVmttLdN2UUOwOYmvo5Afif1O9Op7j7a4PWahERCSxIj74G2OLuze7eASwD5nUrMw/4f57UAHzUzD4+yG0VEZGDECTQjwe2Zxy3pM4FLePAY2a23sx6/ObGzOrMrNHMGnft2hWgWZIrmmkjQWgGTvEIEugtyznvR5kT3X0GyeGdy8zspGwv4u4xd6929+qxY8cGaJbkivaClb6UVUZpHz2RhuY21jQ0KgdOgQsy66YFyExROAHYEbSMu3f+bjWzFSSHgtYcbINlaGgvWOmNZuAUlyA9+nXAVDOrMLMSYAGwsluZlcCFqdk3s4A33f1VMyszs3IAMysD5gIbBrH9IiLShz579O6+z8wuBx4FRgJ3u/tGM1uUur4UWAWcCWwB3gEuTlU/ElhhZp2vFXf3Pw76u5ABUZ55GajO8fpOCxcu1GKqAhJowZS7ryIZzDPPLc147MBlWeo1A8cNsI2SY11Wv4JWwEq/lFVGmd7+fPq48/sdBfrCoZWxAmj1qxy88kgt9UvuSB9n9uylMCjQi8iATV78SPrxzuY2Zk0Zk8fWSHfKXjlMaa685JLm2BcW9eiHKWWllFzJHLPXeH1hUKAfxjQuL7lQHqllG8k59u3NmmNfCDR0IyIScgr0w4jG5SUf1jQ0UjppOqWTpmu8Pk8U6IcR5bCRoVZWGU2vz+ho3aqcOHmiMfphRjlsZCgpJ05hUKAPMaU2kEKjVAn5oUAfYkptIIUkc3NxQBuMDyEF+pDJ7MUnEglNoZSCkTmMA8mhnMwevnr3uaNAHzJaCCXFQgurho4CfQjEYjGuvDmZVKqjdat68VIUui+sUu8+dxToi1BmYAd4d3tyL5cPTZxGybgK9eKl6GSO37+7fQOrV6/uMpFAgX9gLJlKvrBUV1d7Y2NjvptRUDLH3levXg0kA3unsspol/FPkWK1O/FH3t5Unz7u7MicfPLJgIJ+T8xsvbtXZ72mQF+4egvuCuwyXGQGfgX9ng040JtZLXAbya0E73L3Jd2uW+r6mSS3Evyau/81SN1shlOg7z7XPZOCu0hXvQX97obbh8CAAr2ZjQReAk4HWkhuFn6+u2/KKHMmcAXJQH8CcJu7nxCkbjbFGOh7C9i9yTYMk0nBXSS77kM8mfr6EOhNsX5ADDTQzwb+r7t/MXX8XQB3/1FGmV8A9e7+29TxZiAKTO6rbjbl5eV+/PHHB3lvBaOvgN0bBXORwdXbh0B3mbthdf5/fDAfEPm2evXqHgN9kFk344HtGcctJHvtfZUZH7AuAGZWB3R+jO5ZvXr15gBty6UjgNf6W6mzJ9HfOq8/+vN+1xtCB3UvQkr34gOhuBert2c5lwr4/VAI9+ITPV0IEugty7nufwb0VCZI3eRJ9xhQMDlMzayxp0/H4Ub34gO6Fx/QvfhAod+LIIG+BcjMhDUB2BGwTEmAuiIikkNB8tGvA6aaWYWZlQALgJXdyqwELrSkWcCb7v5qwLoiIpJDffbo3X2fmV0OPEpyiuTd7r7RzBalri8FVpGccbOF5PTKi3urm5N3MvgKZhipAOhefED34gO6Fx8o6HtRkAumRERk8GgrQRGRkFOgFxEJOQX6AMzsajNzMzsi323JFzP7sZm9aGZNZrbCzD6a7zYNJTOrNbPNZrbFzIbt5qdmNtHMnjSzF8xso5l9O99tyjczG2lmfzOzh/Pdlp4o0PfBzCaSTOHwr3y3Jc/+BExz9yqSaS2+m+f2DJlUKo87gTOASuB8M6vMb6vyZh/wHXf/DDALuGwY34tO3wZeyHcjeqNA37efAv+HHhZ6DRfu/pi770sdNpBcEzFc1ABb3L3Z3TuAZcC8PLcpL9z91c6Ehe6+m2SAG5/fVuWPmU0AvgTcle+29EaBvhdmdg7wirv/Pd9tKTBfB/6Q70YMoZ5SfAxrZjYZ+Czwlzw3JZ9+RrIjuD/P7ejVsN9hysweB47Kcul7wHXA3KFtUf70di/c/fepMt8j+ef7fUPZtjwLnMpjuDCzjwAPAFe5+1v5bk8+mNlZQKu7rzezaJ6b06thH+jd/bRs581sOlAB/D2Zbp8JwF/NrMbddw5hE4dMT/eik5ldBJwFnOrDawFGkDQgw4aZjSIZ5O9z9wfz3Z48OhE4J5WmvRQYbWa/cfev5LldB9CCqYDMbBtQ7e75zlCXF6kNZH4CnOzuu/LdnqFkZoeQ/AL6VOAVkqk9FhbRKu9Bk9pk6F7gdXe/Ks/NKRipHv3V7n5WnpuSlcboJaifA+XAn8wsYWZL892goZL6ErozlccLwO+GY5BPORH4KvCF1H8HiVSPVgqYevQiIiGnHr2ISMgp0IuIhJwCvYhIyCnQi4iEnAK9iEjIKdCLiIScAr2ISMj9f1s8PMUwwTADAAAAAElFTkSuQmCC\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAgp0lEQVR4nO3df3CU1dn/8fdFgEbSYC2CdvghocXSKHGLIfLD6voDGq3K43QcEFutPprSVq0z1W+pj3W+tU7FqdNWrc+DO5bqtzWFZ1Qso1St1UD9ESG02wgolgYsEWkwVgUlxsj1/WM3YRM2yR3IZjf3fl4zDHvf9zmbszd67cm5z7mOuTsiIhJeQ7LdABERySwFehGRkFOgFxEJOQV6EZGQU6AXEQm5odluQDpHH320T5w4MdvNEBEZNDZs2PCWu49Ody0nA/3EiROpq6vLdjNERAYNM3u9u2sauhERCTkFehGRkFOgFxEJuZwcoxeRcPjoo49obGykpaUl200JjcLCQsaNG8ewYcMC11GgF5GMaWxspLi4mIkTJ2Jm2W7OoOfuNDc309jYSElJSeB6GroRkYxpaWlh1KhRCvL9xMwYNWpUn39DUqAXkYxSkO9fh3I/FehFREJOgV5EpB/F43FWr17d53rRaDRjC0X1MFakDyYufrzj9fYlX8liSyRXxeNx6urqOPfcc7PdlA7q0YtIqP32t7+loqKCSCTCN7/5TV566SXKyspoaWnh/fff54QTTmDjxo3U1NRw2mmnceGFF1JaWsqiRYvYv38/AE899RQzZ85k2rRpXHTRRezduxeA9evXM2vWLE466SQqKip49913ufnmm1mxYgWRSIQVK1bw/vvvc8UVVzB9+nS++MUv8vvf/x6Affv2sWDBAsrKypg/fz779u3L2D1Qj15EBsR1111HPB7v1/eMRCL84he/6Pb6K6+8wooVK3j++ecZNmwY3/72t9myZQsXXHABN910E/v27eNrX/saJ554IjU1Naxbt47Nmzdz3HHHUVlZySOPPEI0GuXWW2/l6aefpqioiNtvv52f/exnLF68mPnz57NixQqmT5/Oe++9x4gRI7jllluoq6vjl7/8JQA33ngjZ555JsuWLeOdd96hoqKCs88+m3vvvZcRI0ZQX19PfX0906ZN69d7kypQoDezSuBOoAC4z92XdFNuOlALzHf3h/pSV0Skv/3pT39iw4YNTJ8+HUj0oseMGcPNN9/M9OnTKSws5K677uooX1FRwaRJkwC4+OKLee655ygsLGTz5s3Mnj0bgNbWVmbOnMmWLVv4zGc+0/HeI0eOTNuGp556ilWrVnHHHXcAiSmn//znP1m7di3XXnstAGVlZZSVlWXmJhAg0JtZAXAPMAdoBNab2Sp335ym3O3Ak32tKyLh11PPO1Pcncsuu4zbbrut0/ldu3axd+9ePvroI1paWigqKgIOnrpoZrg7c+bM4Xe/+12na/X19YGmOro7Dz/8MJ///OcPujZQU0+DjNFXAFvdvcHdW4HlwLw05a4BHgaaDqGuSM6auPjxjj8yuJx11lk89NBDNDUlwtLbb7/N66+/TlVVFT/+8Y+55JJL+P73v99Rft26dWzbto39+/ezYsUKTj31VGbMmMHzzz/P1q1bAfjggw947bXXmDJlCjt37mT9+vUA7Nmzh7a2NoqLi9mzZ0/He375y1/m7rvvxt0B+Otf/wrAaaedxoMPPgjAxo0bqa+vz9h9CDJ0MxbYkXLcCJySWsDMxgIXAmcC0/tSN+U9qoAqgAkTJgRolkjmxWIxdlXfnfbaqPgTFEcqO441Cyf3lJaWcuuttzJ37lz279/PsGHDmDdvHkOHDmXhwoV8/PHHzJo1i2eeeYYhQ4Ywc+ZMFi9ezMsvv9zxYHbIkCHcf//9XHzxxXz44YcA3HrrrRx//PGsWLGCa665hn379nHEEUfw9NNPc8YZZ7BkyRIikQg/+MEP+OEPf8h1111HWVkZ7s7EiRN57LHH+Na3vsXll19OWVkZkUiEioqKjN2HIIE+3e8W3uX4F8D33f3jLr+KBKmbOOkeA2IA5eXlacuIDLTq6mpam7YxfEznvCKtTdsAOgV6yU3z589n/vz5aa8VFBTw0ksvAVBTU8OIESNYsWLFQeXOPPPMjp57qunTp1NbW3vQ+a5l77333oPKHHHEESxfvjzQZzhcQQJ9IzA+5XgcsLNLmXJgeTLIHw2ca2ZtAeuK5IxYLMa1tx7owbcH+WMXdp5DsKt6Ma1N29hVvfhA3UlvUFVVNWBtFQkqyBj9emCymZWY2XBgAbAqtYC7l7j7RHefCDwEfNvdHw1SVySXtPfg2w0fU0JRafSgckWl0U69/NambVRXVw9EEyVDotEojz32WLabkRG99ujdvc3MriYxm6YAWObum8xsUfL60r7W7Z+mi2RGuh58V8WRyk7DNqk9e5FcE2gevbuvBlZ3OZc2wLv7N3qrKyIiA0cpECTvxWIxotEo0Wj0sFZuxuPxjveJxWL910CRw6RAL3mvurq6I8BHIpG0Y/K9KSqNEolEgETA13i95BLluhEhEeBramoADmlhVHGkku0kxuxbGjRe353+XnQWZO3C9u3bOe+889i4cWOn8zfffDOnnXYaZ599dtp6jz76KMcffzylpaX90tZsUqAXAWobmrXyNc/ccsstPV5/9NFHOe+88/ol0Le1tTF0aPbCrYZuRCT0Pv74Y6666ipOOOEE5s6dy759+/jGN77BQw89BMDixYspLS2lrKyM66+/nhdeeIFVq1Zxww03EIlE+Mc//kE8HmfGjBmUlZVx4YUX8u9//xtILI4qKytj5syZ3HDDDZx44okA3H///Vx00UWcf/75zJ07l71793LWWWcxbdo0pk6d2pGuePv27UyZMoUrr7ySE088kUsuuYSnn36a2bNnM3nyZNatW3fYn1+BXkRC7+9//zvf+c532LRpE5/61Kd4+OGHO669/fbbrFy5kk2bNlFfX89NN93ErFmzuOCCC/jpT39KPB7ns5/9LJdeeim333479fX1TJ06lR/96EcAXH755SxdupQXX3yRgoKCTj/3xRdf5IEHHuCZZ56hsLCQlStX8pe//IVnn32W733vex35b7Zu3cp3v/td6uvrefXVV6murua5557jjjvu4Cc/+clhf34FeslL/TXTpjuagZNbSkpKOh6Wn3zyyWzfvr3j2siRIyksLOTKK6/kkUceYcSIEQfVf/fdd3nnnXc4/fTTAbjssstYu3Yt77zzDnv27GHWrFkALFy4sFO9OXPm8OlPfxpIZLG88cYbKSsr4+yzz+aNN97gX//6V0f7pk6dypAhQzjhhBM466yzMDOmTp3aqa2HSoFe8lJ1dTVra+uobWimZeT4Q5pp0x3NwMk9n/jEJzpeFxQU0NbW1nE8dOhQ1q1bx1e/+lUeffRRKiuD5y9q75F3pz39McCDDz7I7t272bBhA/F4nGOOOYaWlpaD2jdkyJCO4yFDhnRq66HSw1jJW0FWwB6K4kglNUsS+XKi0Wi/v7/0r7179/LBBx9w7rnnMmPGDD73uc8BdEo3fOSRR3LUUUfx5z//mS996Uv85je/4fTTT+eoo46iuLiY2tpaZsyY0WOSsnfffZcxY8YwbNgwnn32WV5//fUB+XygQC8iAygXUznv2bOHefPm0dLSgrvz85//HIAFCxZw1VVXcdddd/HQQw/xwAMPsGjRIj744AMmTZrEr3/9awB+9atfcdVVV1FUVEQ0GuXII49M+3MuueQSzj//fMrLy4lEIkyZMmXAPqP19qtHNpSXl3tdXV22myEhFo1GqW1ozkiPHg4EtPYeffsc/Xzzyiuv8IUvfCHbzciovXv38slPfhKAJUuW8Oabb3LnnXdm9Gemu69mtsHdy9OVV49e8kbqPPldDc0D8rN2NTQzY9KojP4sya7HH3+c2267jba2No477jjuv//+bDfpIAr0IiKHoaeNTXKFAr3kha5bAqbbNSpT2qdatlu4cGFebVDi7gO2CXY+OJThdk2vlLwQdEOR/lZUGqVl5HhqG5qpbWhmbW1dXk23LCwspLm5+ZCCkxzM3WlubqawsLBP9dSjl7yRqemUPcn3DUrGjRtHY2Mju3fvznZTQqOwsJBx48b1qU6gQG9mlcCdJHaJus/dl3S5Pg/4MbAfaAOuc/fnkte2A3uAj4G27p4Ki0j4DBs2jJKSgRkik+71GujNrAC4B5hDYrPv9Wa2yt03pxT7E7DK3d3MyoD/BVIniZ7h7m/1Y7tFRCSgIGP0FcBWd29w91ZgOTAvtYC77/UDg3BFgAbkRERyRJBAPxbYkXLcmDzXiZldaGavAo8DV6RccuApM9tgZt1ONTCzKjOrM7M6jeeJiPSfIIE+3byog3rs7r7S3acA/0FivL7dbHefBpwDfMfMTkv3Q9w95u7l7l4+evToAM0SEZEgggT6RmB8yvE4YGd3hd19LfBZMzs6ebwz+XcTsJLEUJCIiAyQIIF+PTDZzErMbDiwAFiVWsDMPmfJFRFmNg0YDjSbWZGZFSfPFwFzgc4bN4pkSKZzzh8q5aqXgdbrrBt3bzOzq4EnSUyvXObum8xsUfL6UuCrwKVm9hGwD5ifnIFzDLAy+R0wFKh29ycy9FlEOqmuriYejxOJRIhEIrxcODXbTaKoNMrUlpcBOr588mmVrGSHsldKaHXNHJkrm38rs6VkgrJXSt6qbWjOmQAvki0K9CIDTCmMZaApqZmISMgp0IuIhJwCvYhIyCnQi2RR+8NiPTCWTFKgFxEJOc26kVCJxWIdOzjF43EYOb7nClnW2rStYzOS2KQ3tHhKMkI9egmV9tWwAJFIZEC2CzxURaXRjn1rW5u25dUWgzKw1KOX0IlEIjm3Gjad1G0G822LQRlYCvQSOloNK9KZhm5EREJOgV5EJOQU6EVEQk6BXkQk5BToRURCLlCgN7NKM9tiZlvN7KB5YGY2z8zqzSxuZnVmdmrQuiIiklm9BnozKwDuAc4BSoGLzay0S7E/ASe5ewS4ArivD3VFRCSDgsyjrwC2unsDgJktB+YBm9sLuPvelPJFgAetK3I4UlMewOBIe9CdtbV1FE44sK/tXTddo5QI0i+CDN2MBXakHDcmz3ViZhea2avA4yR69YHrJutXJYd96nbv3h2k7SKdUh5A7qc96E5qOgRQSgTpX0F69Jbm3EE7irv7SmClmZ0G/Bg4O2jdZP0YEIPE5uAB2iUCQMvI8WyfcUPHcXEW23KoUtMhgFIiSP8K0qNvBFJ/Fx4H7OyusLuvBT5rZkf3ta6IiPS/IIF+PTDZzErMbDiwAFiVWsDMPmdmlnw9DRgONAepKyIimdXr0I27t5nZ1cCTQAGwzN03mdmi5PWlwFeBS83sI2AfMN/dHUhbN0OfRURE0giUvdLdVwOru5xbmvL6duD2oHVFRGTgaGWsiEjIKdCLiIScAr2ISMgp0IuIhJwCvYhIyCnQy6ATi8WIRqNEo9FO6Q/CJh6Pd3zOWCyW7ebIIKZAL4NOan6bwZrbpjdFpVEikQiQCPjKeyOHI9A8epFcE4lEqKmpAWDi4sez25gMKI5UUrPkbgCi0Wh2GyODngK9DEq1Dc2hDPAimaBAL5Kj2r/IdjU0M2PSqCy3RgYzjdGLiIScAr2ISMgp0IuIhJwCvYhIyCnQi4iEnAK9iEjIBZpeaWaVwJ0kdom6z92XdLl+CfD95OFe4Fvu/rfkte3AHuBjoM3dy/un6ZJPUufM72pozmJLRAafXgO9mRUA9wBzSGz2vd7MVrn75pRi24DT3f3fZnYOEANOSbl+hru/1Y/tljwSi8XYVX13x3Fr0zaGjynJYosGXnveG4CFCxdSVVWV3QbJoBJk6KYC2OruDe7eCiwH5qUWcPcX3P3fycNaYFz/NlPyWXV1Na1N2zqOh48pCWV+m+4o740criBDN2OBHSnHjXTurXf1n8AfUo4deMrMHLjX3ZWGT/ps+JgSjl24pPeCIVQcqWQ7lQC0NCzOcmtkMAoS6C3NOU9b0OwMEoH+1JTTs919p5mNAf5oZq+6+9o0dauAKoAJEyYEaJaIiAQRZOimERifcjwO2Nm1kJmVAfcB89y942mZu+9M/t0ErCQxFHQQd4+5e7m7l48ePTr4JxARkR4FCfTrgclmVmJmw4EFwKrUAmY2AXgE+Lq7v5ZyvsjMittfA3OBjf3VeBER6V2vQzfu3mZmVwNPkpheuczdN5nZouT1pcDNwCjgv80MDkyjPAZYmTw3FKh29ycy8klERCStQPPo3X01sLrLuaUpr68ErkxTrwE46TDbKCIih0ErY0VEQk6BXkQk5LTDlOSs1B2WROTQqUcvIhJy6tFLTkrNb5OPuW16kpr3BpT7RnqnHr3kpNT8NvmW26YnqXlvQLlvJBj16CVn5XN+m+6k5r0B5b6RYNSjFxEJOQV6EZGQU6AXEQk5BXoRkZBToBcRCTkFehGRkFOgFxEJOQV6EZGQU6AXEQm5QIHezCrNbIuZbTWzg5bimdklZlaf/POCmZ0UtK5Iu1gsRjQaJRqNEo/Hs92cQaM99000GiUWi2W7OZKDeg30ZlYA3AOcA5QCF5tZaZdi24DT3b0M+DEQ60NdESCR36Y9wEciEeW3CSA1943y3kh3guS6qQC2JrcFxMyWA/OAze0F3P2FlPK1wLigdUVStYwcz/YZNwBQnOW2DAapuW+U90a6E2ToZiywI+W4MXmuO/8J/OEQ64qISD8L0qO3NOc8bUGzM0gE+lMPoW4VUAUwYcKEAM0SEZEggvToG4HxKcfjgJ1dC5lZGXAfMM/dm/tSF8DdY+5e7u7lo0ePDtJ2EREJIEigXw9MNrMSMxsOLABWpRYwswnAI8DX3f21vtQVEZHM6nXoxt3bzOxq4EmgAFjm7pvMbFHy+lLgZmAU8N9mBtCW7J2nrZuhzyIiImkE2mHK3VcDq7ucW5ry+krgyqB1RURk4GhlrIhIyCnQi4iEnDYHl6yJxWKdVnLG43EYOb77CtKr9nQIAAsXLqSqqiq7DZKcoB69ZE1qygNQ2oPDVVQapWXkeGobmllbW6d0CNJBPXrJqkgkQk1NTcfxxMWPZ68xg1xxpJLiSCIdwq5qpUOQA9SjFxEJOfXoJatqG5rVixfJMPXoRURCToFeRCTkFOhFREJOgV5EJOQU6EVEQk6BXkQk5DS9UgZUatoDpTzIrNR0CKCUCPlMgV4GVHV1NWtr6xg+pgRGjlfKgwwpKo0yteXljuP2VBMK9PlJgV4G3PAxJRy7cEm2mxFqxZFKapbc3XGc2rOX/BNojN7MKs1si5ltNbODkmiY2RQze9HMPjSz67tc225mL5tZ3Mzq+qvhIiISTK89ejMrAO4B5pDY7Hu9ma1y980pxd4GrgX+o5u3OcPd3zrMtopIH6SmltjV0MyMSaOy2BrJpiA9+gpgq7s3uHsrsByYl1rA3ZvcfT3wUQbaKCIihyFIoB8L7Eg5bkyeC8qBp8xsg5npSZCIyAAL8jDW0pzzPvyM2e6+08zGAH80s1fdfe1BPyTxJVAFMGHChD68vYiI9CRIj74RSJ3sPA7YGfQHuPvO5N9NwEoSQ0HpysXcvdzdy0ePHh307UVEpBdBAv16YLKZlZjZcGABsCrIm5tZkZkVt78G5gIbD7WxIiLSd70O3bh7m5ldDTwJFADL3H2TmS1KXl9qZscCdcBIYL+ZXQeUAkcDK82s/WdVu/sTGfkkIiKSVqAFU+6+Gljd5dzSlNe7SAzpdPUecNLhNFAGv1gsxrW3JhbvtDZtS6yKlQG3traOwglTAbjrpmu0SjaPKKmZZFx1dTWtTduAxKpYpT0YeEWl0Y4v2NambR35hiQ/KAWCDAilPciu4kglxZFKAHZVH7S4XUJOPXoRkZBToBcRCTkFehGRkFOgFxEJOQV6EZGQU6AXEQk5Ta+UjOiaC11yS21Dc6d/o+1LvpLF1kimKdBLv4vFYuyqPrCNnVbD5p7Wpm2d5tPHJr2hlbIhpqEb6XepK2FBq2FzTeoqWdBK2XygHr1khFbC5q7UVbKglbL5QD16EZGQU6AXEQk5BXoRkZBToBcRCTkFehGRkAsU6M2s0sy2mNlWMzvoEb2ZTTGzF83sQzO7vi91RUQks3oN9GZWANwDnENiH9iLzay0S7G3gWuBOw6hroRALBYjGo0SjUaJx+PZbo70UTwe7/j3i8Vi2W6O9LMgPfoKYKu7N7h7K7AcmJdawN2b3H098FFf60o4VFdXs7a2jtqGZlpGjtcCqUGkqDRKJBIBEgFfi6fCJ8iCqbHAjpTjRuCUgO8fuK6ZVQFVABMmTAj49pJLtEhqcCqOVLKdxAKqlgaNroZRkB69pTnnAd8/cF13j7l7ubuXjx49OuDbi4hIb4IE+kZgfMrxOGBnwPc/nLoiItIPggT69cBkMysxs+HAAmBVwPc/nLoiItIPeh2jd/c2M7saeBIoAJa5+yYzW5S8vtTMjgXqgJHAfjO7Dih19/fS1c3QZxERkTQCZa9099XA6i7nlqa83kViWCZQXRHJXambkmhDknBQmmI5ZLFYrGMqXjweh5Hje64gg0LqpiTakCQclAJBDpnmzodP6qYk2pAkPNSjl8OiufPhkropiTYkCQ/16EVEQk6BXkQk5BToRURCToFeRCTkFOglsFgsRuGEqR1/1tbWZbtJkmFra+s6/ZsrhfHgpEAvgVVXV9PatK3jePiYEk2pDLHUqZag6ZaDmaZXSp9oOmX+SJ1qCZpuOZipRy8iEnIK9CIiIadALyIScgr00iNt+i2ptIn44KSHsdKj9sRlw8eUgBKX5bWi0ihTW14G6PjSV2bLwUGBXnqlmTYC2kR8MAs0dGNmlWa2xcy2mtlB/8KWcFfyer2ZTUu5tt3MXjazuJlphY2IyADrtUdvZgXAPcAcEpt9rzezVe6+OaXYOcDk5J9TgP9J/t3uDHd/q99aLSIigQXp0VcAW929wd1bgeXAvC5l5gH/zxNqgU+Z2Wf6ua0yQPQAVoLQg9nBI0igHwvsSDluTJ4LWsaBp8xsg5l1++TGzKrMrM7M6nbv3h2gWZIp2jlKelNUGqVl5HhqG5pZW1un1Ag5LsjDWEtzzvtQZra77zSzMcAfzexVd197UGH3GBADKC8v7/r+MsD0AFZ6op2oBpcgPfpGIHXX53HAzqBl3L397yZgJYmhIBERGSBBevTrgclmVgK8ASwAFnYpswq42syWk3gI+667v2lmRcAQd9+TfD0XuKX/mi8iuaC2oZmJix/vON6+5CtZbI101Wugd/c2M7saeBIoAJa5+yYzW5S8vhRYDZwLbAU+AC5PVj8GWGlm7T+r2t2f6PdPIYclFot1GmONx+Mwcnz3FUS6aG3a1mkIJzbpDS2myiHmnnvD4eXl5V5Xpyn3AyUajR5Y/ZpUVBrtlKJWpDt74k/w/uaajuPWpm2cNqOcmpqabutI/zOzDe5enu6aVsYKoIevcuiUtz73KamZiEjIKdDnKS2KkkzSYqrcoqGbPKWslJIpRaVR3t9cQ21Dc8cew3owm10K9HlM4/KSCVpMlXs0dJNHNFwj2bC2to7CCVMpnDBVwzhZokCfR6qrqzsCfCQS0XCNZFxRabRj2m5r0zblxMkSDd3kmZaR49k+4wYAirPcFgk/DePkBgX6ENOKV8k17bNx2i1cuFAPageAAn2IdZpZA5pdI1mVOhsH0IycAaRAHzKpvfh4PK6ZNZIz0q2gTe3hq3efOQr0IaP58TJYaL79wFGgDwH14mUw6vqgVr37zFGgH4S6PmRds2YNAJ8Yf6J68TIopfbuP9yxkTVr1nT6b1yB//AoTfEgkRrcOwX2JKUVlrDomvb4wx0bATj99NMBBf3u9JSmWIE+h/UU3BXYJV+kBn4F/e4ddqA3s0rgThI7TN3n7ku6XLfk9XNJ7DD1DXf/S5C66eRToO86DJOqPbi3/0f9cuFUBXfJaz0F/a7y7UvgsAK9mRUArwFzSGwCvh642N03p5Q5F7iGRKA/BbjT3U8JUjedwRjoewrYPUk3DJNKPXeR9LoO8QDMmDQKOLiT1BeD9QvicAP9TOD/uvuXk8c/AHD321LK3AvUuPvvksdbgCgwsbe66RQXF/vJJ58c5LNlXfvij3S9i/ZrvVEwF+lf6b4Egmj//7i949X+xTEYrFmz5rC2EhwL7Eg5biTRa++tzNiAdQEwsyqg/Wt075o1a7YEaFsmHQ281ddK7T2Jvvhwx0befvKXfa43gA7pXoSU7sUBob0X7QF/zY5eCh6QC/fiuO4uBAn0luZc118DuisTpG7ipHsMyJkcpmZW1923Y77RvThA9+IA3YsDcv1eBAn0jUBqJqxxwM6AZYYHqCsiIhkUJB/9emCymZWY2XBgAbCqS5lVwKWWMAN4193fDFhXREQyqNcevbu3mdnVwJMkpkguc/dNZrYoeX0psJrEjJutJKZXXt5T3Yx8kv6XM8NIOUD34gDdiwN0Lw7I6XuRkwumRESk/2grQRGRkFOgFxEJOQX6AMzsejNzMzs6223JFjP7qZm9amb1ZrbSzD6V7TYNJDOrNLMtZrbVzPJ281MzG29mz5rZK2a2ycy+m+02ZZuZFZjZX83ssWy3pTsK9L0ws/EkUjj8M9ttybI/Aie6exmJtBY/yHJ7Bkwylcc9wDlAKXCxmZVmt1VZ0wZ8z92/AMwAvpPH96Ldd4FXst2InijQ9+7nwP+hm4Ve+cLdn3L3tuRhLYk1EfmiAtjq7g3u3gosB+ZluU1Z4e5vticsdPc9JALc2Oy2KnvMbBzwFeC+bLelJwr0PTCzC4A33P1v2W5LjrkC+EO2GzGAukvxkdfMbCLwReClLDclm35BoiO4P8vt6FHe7zBlZk8Dx6a59F/AjcDcgW1R9vR0L9z998ky/0Xi1/cHB7JtWRY4lUe+MLNPAg8D17n7e9luTzaY2XlAk7tvMLNolpvTo7wP9O5+drrzZjYVKAH+lki3zzjgL2ZW4e67BrCJA6a7e9HOzC4DzgPO8vxagBEkDUjeMLNhJIL8g+7+SLbbk0WzgQuSadoLgZFm9lt3/1qW23UQLZgKyMy2A+Xunu0MdVmR3EDmZ8Dp7r472+0ZSGY2lMQD6LOAN0ik9lg4iFZ595vkJkMPAG+7+3VZbk7OSPbor3f387LclLQ0Ri9B/RIoBv5oZnEzW5rtBg2U5EPo9lQerwD/m49BPmk28HXgzOR/B/Fkj1ZymHr0IiIhpx69iEjIKdCLiIScAr2ISMgp0IuIhJwCvYhIyCnQi4iEnAK9iEjI/X92ak1yD3ykGAAAAABJRU5ErkJggg==\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(-np.log(p**2)))\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
|
|
}
|