931 lines
94 KiB
Plaintext
931 lines
94 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": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD7CAYAAACRxdTpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAZKElEQVR4nO3df3BV9Z3/8ec7IRrBoJYfXb9EJDgoRpJmbUDALlz5ZaQqw3Y7SGhVWsjYKtWZ1q/oUGyVqTA6LeKPRcYibmskO4CWVXZVqkBRKYT2bvjV2BRQI1oiCAVNwMB7/0iEGELuuckNNzl5PWYyzTnnfe99e6a87ifnfu7nmLsjIiLhkpLsBkREJPEU7iIiIaRwFxEJIYW7iEgIKdxFREJI4S4iEkIxw93MFpvZXjPbGqNusJkdM7N/S1x7IiLSEkFG7kuAguYKzCwVmAe8koCeRESklbrEKnD3dWbWL0bZDGA5MDjoC/fs2dP79Yv1tCIi0tDmzZs/dvdesepihnssZtYHmAiMIo5w79evH6Wlpa19eRGRTsXM3g1Sl4gPVOcD97j7sViFZlZkZqVmVlpVVZWAlxYRkaa0euQO5ANLzQygJzDezGrd/cXGhe6+CFgEkJ+fr0VtRETaSKvD3d2zvvjdzJYALzUV7CIicubEDHczex6IAD3NrBK4H0gDcPeFbdqdiLR7n3/+OZWVldTU1CS7lVBJT08nMzOTtLS0Fj0+yGyZyUGfzN1vbVEXItJhVVZWkpGRQb9+/ai/PCut5O7s27ePyspKsrKyYj+gCfqGqoi0Sk1NDT169FCwJ5CZ0aNHj1b9NaRwF5FWU7AnXmvPqcJdRCSEFO4iIi0QjUZZtWpV3I+LRCJn5AuciZjnLq3Ub+bLgWt3z/1mG3YiIkFFo1FKS0sZP358sltpksI9yRYtWsRHxY8Fr+//AUVFRW3YkUjH9Nvf/pYFCxZw9OhRrrrqKr73ve8xffp0Nm7cyLFjxxgyZAglJSV8/PHHzJ49mx49elBeXs6IESN48sknSUlJ4dVXX+X+++/nyJEjXHLJJTzzzDOce+65bNq0iTvvvJNPP/2Us88+m9dee43Zs2dTXV3N+vXruffee7n++uuZMWMGW7Zsoba2lp/97GdMmDCB6upqpk6dyvbt27n88suprq4+I+dD4Z5kxcXFHN27i7N6x57udHTvLoqLixXu0m7dddddRKPRhD5nXl4e8+fPb7Zmx44dlJSU8Oabb5KWlsYPf/hDysvLufHGG5k1axbV1dV85zvfYdCgQaxZs4aNGzeyfft2Lr74YgoKClixYgWRSIQ5c+awevVqunXrxrx58/jlL3/JzJkzmTRpEiUlJQwePJh//OMfdO3alQceeIDS0lIef/xxAO677z5GjRrF4sWLOXDgAEOGDGHMmDE89dRTdO3albKyMsrKyrjyyisTen5OR+HeDpzVO4t/Kpwbs+6j4plnoBuRjuf3v/89mzdvZvDgurULq6ur6d27N7Nnz2bw4MGkp6ezYMGCE/VDhgyhf//+AEyePJn169eTnp7O9u3bufrqqwE4evQow4YNo7y8nAsvvPDEc3fv3r3JHl599VVWrlzJI488AtRNEX3vvfdYt24dP/rRjwDIzc0lNze3bU5CIwp3EUmYWCPstuLu3HLLLTz00ENf2v/RRx9x+PBhPv/8c2pqaujWrRtw6jRDM8PdGTt2LM8///yXjpWVlQWalujuLF++nMsuu+yUY8mYKqrZMiLS4Y0ePZply5axd+9eAPbv38+7775LUVERDz74IFOmTOGee+45Ub9x40Z27drF8ePHKSkp4Rvf+AZDhw7lzTffpKKiAoDPPvuMd955h4EDB7Jnzx42bdoEwKFDh6itrSUjI4NDhw6deM5rr72Wxx57DPe6NRH//Oc/AzBixAiee+45ALZu3UpZWVnbnxA0cheREMjOzmbOnDmMGzeO48ePk5aWxoQJE+jSpQuFhYUcO3aM4cOH8/rrr5OSksKwYcOYOXMmW7ZsYcSIEUycOJGUlBSWLFnC5MmTOXLkCABz5szh0ksvpaSkhBkzZlBdXc0555zD6tWrueaaa5g7dy55eXnce++9/PSnP+Wuu+4iNzcXd6dfv3689NJL/OAHP2Dq1Knk5uaSl5fHkCFDzsg5UbiLSChMmjSJSZMmNXksNTWVP/7xjwCsWbOGrl27UlJSckrdqFGjTozQGxo8eDAbNmw4ZX/j2qeeeuqUmnPOOYelS5cG+m9IJF2WEREJIY3cRaRTiUQiRCKRZLfR5jRyFxEJIYW7iEgIKdxFREJI4S4iEkL6QFVEEiqeVU6DCLIS6u7du7n++uvZunXrl/bPnj2bESNGMGbMmCYf9+KLL3LppZeSnZ2dkF7bE43cRSS0HnjggdMGO9SF+/bt2xPyWrW1tQl5nkRRuItIKBw7dozp06dzxRVXMG7cOKqrq7n11ltZtmwZADNnziQ7O5vc3Fx+8pOf8NZbb7Fy5Uruvvtu8vLy+Nvf/kY0GmXo0KHk5uYyceJEPvnkE6Duy0q5ubkMGzaMu+++m0GDBgGwZMkSvv3tb3PDDTcwbtw4Dh8+zOjRo7nyyivJycnhd7/7HVD3l8XAgQOZNm0agwYNYsqUKaxevZqrr76aAQMGsHHjxoSfD4W7iITCX//6V26//Xa2bdvG+eefz/Lly08c279/Py+88ALbtm2jrKyMWbNmMXz4cG688UYefvhhotEol1xyCTfffDPz5s2jrKyMnJwcfv7znwMwdepUFi5cyNtvv01qauqXXvftt9/m2Wef5fXXXyc9PZ0XXniBP/3pT7zxxhv8+Mc/PrHWTEVFBXfeeSdlZWX85S9/obi4mPXr1/PII4/wi1/8IuHnI2a4m9liM9trZltPc3yKmZXV/7xlZl9LeJciIjFkZWWRl5cHwNe//nV279594lj37t1JT09n2rRprFixgq5du57y+IMHD3LgwAFGjhwJwC233MK6des4cOAAhw4dYvjw4QAUFhZ+6XFjx47lK1/5ClC3MuR9991Hbm4uY8aM4YMPPuDvf//7if5ycnJISUnhiiuuYPTo0ZgZOTk5X+o1UYKM3JcABc0c3wWMdPdc4EFgUQL6EhGJy9lnn33i99TU1C9dA+/SpQsbN27kW9/6Fi+++CIFBc1F2pd9MfI+nS+WEQZ47rnnqKqqYvPmzUSjUb761a9SU1NzSn8pKSkntlNSUtrken3McHf3dcD+Zo6/5e6f1G9uADIT1JuISEIcPnyYgwcPMn78eObPn3/iblENl+0977zzuOCCC/jDH/4AwG9+8xtGjhzJBRdcQEZGxomFw5pbBOzgwYP07t2btLQ03njjDd599922/Q9rRqKnQn4f+O/THTSzIqAIoG/fvgl+aRFpD9rjTdwPHTrEhAkTqKmpwd351a9+BcBNN93E9OnTWbBgAcuWLePZZ5/ltttu47PPPqN///4888wzAPz6179m+vTpdOvWjUgkwnnnndfk60yZMoUbbriB/Px88vLyGDhw4Bn7b2zMYv3JAWBm/YCX3H1QMzXXAE8C33D3fbGeMz8/30tLS+NoNZwikQgbdu4LfJu9of17sGbNmrZvTCSgHTt2cPnllye7jTZ1+PBhzj33XADmzp3Lhx9+yKOPPtrmr9vUuTWzze6eH+uxCRm5m1ku8DRwXZBgFxHpSF5++WUeeughamtrufjii1myZEmyW4qp1eFuZn2BFcB33f2d1rckItK+NHcjkPYqZrib2fNABOhpZpXA/UAagLsvBGYDPYAn628CWxvkTwYRCQ93T8pNoMMsyCXz5sQMd3efHOP4NGBaq7oQkQ4rPT2dffv20aNHDwV8grg7+/btIz09vcXPoYXDRKRVMjMzqayspKqqKtmthEp6ejqZmS2fWa5wF5FWSUtLIysrK9ltSCNaW0ZEJIQU7iIiIaRwFxEJIYW7iEgIKdxFREJI4S4iEkIKdxGREFK4i4iEkMJdRCSEFO4iIiGkcBcRCSGFu4hICCncRURCSOEuIhJCCncRkRBSuIuIhJDCXUQkhHQnpg4mGo0SiUQC1RYWFlJUVNS2DYlIu6Rw70C6ZUfIqdkSqDYajQIo3EU6KYV7B5KRV8BuCgLV1uyc2cbdiEh7FvOau5ktNrO9Zrb1NMfNzBaYWYWZlZnZlYlvU0RE4hHkA9Ul0Oxw8TpgQP1PEfDvrW9LRERaI2a4u/s6YH8zJROA//A6G4DzzezCRDUoIiLxS8RUyD7A+w22K+v3ncLMisys1MxKq6qqEvDSIiLSlESEuzWxz5sqdPdF7p7v7vm9evVKwEuLiEhTEhHulcBFDbYzgT0JeF4REWmhRIT7SuDm+lkzQ4GD7v5hAp5XRERaKOY8dzN7HogAPc2sErgfSANw94XAKmA8UAF8Bkxtq2ZFRCSYmOHu7pNjHHfg9oR1JCIiraaFw0REQkjhLiISQgp3EZEQUriLiISQwl1EJIQU7iIiIaRwFxEJIYW7iEgIKdxFREJI4S4iEkIKdxGREFK4i4iEkMJdRCSEFO4iIiGkcBcRCSGFu4hICCncRURCSOEuIhJCCncRkRBSuIuIhJDCXUQkhBTuIiIhFCjczazAzMrNrMLMZjZx/Dwz+y8z+18z22ZmUxPfqoiIBNUlVoGZpQJPAGOBSmCTma109+0Nym4Htrv7DWbWCyg3s+fc/WibdC2BrNtQSnrfnMD1C2bNoKioqA07EpEzJcjIfQhQ4e4768N6KTChUY0DGWZmwLnAfqA2oZ1KXLplRzird1bg+qN7d1FcXNyGHYnImRRz5A70Ad5vsF0JXNWo5nFgJbAHyAAmufvxxk9kZkVAEUDfvn1b0q8ElJFXQEZeQeD6j4pPudomIh1YkJG7NbHPG21fC0SB/wfkAY+bWfdTHuS+yN3z3T2/V69ecbYqIiJBBQn3SuCiBtuZ1I3QG5oKrPA6FcAuYGBiWhQRkXgFCfdNwAAzyzKzs4CbqLsE09B7wGgAM/sqcBmwM5GNiohIcDGvubt7rZndAbwCpAKL3X2bmd1Wf3wh8CCwxMy2UHcZ5x53/7gN+xYRkWYE+UAVd18FrGq0b2GD3/cA4xLbmoiItJS+oSoiEkIKdxGREFK4i4iEkMJdRCSEFO4iIiGkcBcRCSGFu4hICCncRURCSOEuIhJCCncRkRBSuIuIhJDCXUQkhBTuIiIhpHAXEQkhhbuISAgp3EVEQkjhLiISQgp3EZEQUriLiISQwl1EJIQC3SBbOodoNEokEglUW1hYSFFRUds2JCItpnAXALplR8ip2RKoNhqNAijcRdqxQOFuZgXAo0Aq8LS7z22iJgLMB9KAj919ZMK6lDaXkVfAbgoC1dbsnNnG3YhIa8UMdzNLBZ4AxgKVwCYzW+nu2xvUnA88CRS4+3tm1ruN+hURkQCCfKA6BKhw953ufhRYCkxoVFMIrHD39wDcfW9i2xQRkXgECfc+wPsNtivr9zV0KXCBma0xs81mdnNTT2RmRWZWamalVVVVLetYRERiChLu1sQ+b7TdBfg68E3gWuCnZnbpKQ9yX+Tu+e6e36tXr7ibFRGRYIJ8oFoJXNRgOxPY00TNx+7+KfCpma0Dvga8k5AuRUQkLkFG7puAAWaWZWZnATcBKxvV/A74FzPrYmZdgauAHYltVUREgoo5cnf3WjO7A3iFuqmQi919m5ndVn98obvvMLP/AcqA49RNl9zalo2LiMjpBZrn7u6rgFWN9i1stP0w8HDiWhMRkZbS2jIiIiGkcBcRCSGFu4hICCncRURCSOEuIhJCCncRkRBSuIuIhJDCXUQkhBTuIiIhpHAXEQkh3UNVWmTdhlLS++YEql0wa4butypyhmnkLnHrlh3hrN5ZgWqP7t1FcXFxG3ckIo1p5C5xy8grICMv2M20PyrWzbRFkkEjdxGREFK4i4iEkMJdRCSEFO4iIiGkcBcRCSGFu4hICCncRURCSOEuIhJCCncRkRAKFO5mVmBm5WZWYWan/cqhmQ02s2Nm9m+Ja1FEROIVM9zNLBV4ArgOyAYmm1n2aermAa8kukkREYlPkJH7EKDC3Xe6+1FgKTChiboZwHJgbwL7ExGRFggS7n2A9xtsV9bvO8HM+gATgYXNPZGZFZlZqZmVVlVVxduriIgEFCTcrYl93mh7PnCPux9r7oncfZG757t7fq9evQK2KCIi8Qqy5G8lcFGD7UxgT6OafGCpmQH0BMabWa27v5iIJqVji0ajRCKRQLWFhYW6sYdIAgQJ903AADPLAj4AbgIKGxa4+4k7N5jZEuAlBbtA3Y09cmq2BKqNRqMACneRBIgZ7u5ea2Z3UDcLJhVY7O7bzOy2+uPNXmeXzi0jr4DdBLuxR81O3dhDJFEC3YnJ3VcBqxrtazLU3f3W1rclIiKtoW+oioiEkMJdRCSEFO4iIiGkcBcRCSGFu4hICCncRURCSOEuIhJCCncRkRBSuIuIhFCgb6iKnClaZEwkMRTu0m50y47w6fY1bNi5L2bt0b27AC0yJnI6CndpNzLyCsjIC7bI2EfFWmRMpDm65i4iEkIKdxGREFK4i4iEkMJdRCSEFO4iIiGkcBcRCSGFu4hICGmeu3RY+jaryOkp3KVD6pYdIadmS6DaaDQK6Nus0rko3KVDysgrYDfBvs1as1PfZpXOJ9A1dzMrMLNyM6sws1P+pZjZFDMrq/95y8y+lvhWRUQkqJjhbmapwBPAdUA2MNnMshuV7QJGunsu8CCwKNGNiohIcEFG7kOACnff6e5HgaXAhIYF7v6Wu39Sv7kByExsmyIiEo8g4d4HeL/BdmX9vtP5PvDfTR0wsyIzKzWz0qqqquBdiohIXIKEuzWxz5ssNLuGunC/p6nj7r7I3fPdPb9Xr17BuxQRkbgEmS1TCVzUYDsT2NO4yMxygaeB69w99t0WRM6geObEg+bFS8cXJNw3AQPMLAv4ALgJKGxYYGZ9gRXAd939nYR3KdIK8dzhCXSXJwmHmOHu7rVmdgfwCpAKLHb3bWZ2W/3xhcBsoAfwpJkB1Lp7ftu1LRJcPHd4At3lScIh0JeY3H0VsKrRvoUNfp8GTEtsayIi0lJaOExEJIQU7iIiIaS1ZUSaoBUnpaNTuIs0Es/sGs2skfZK4S7SSDyzazSzRtorXXMXEQkhhbuISAjpsoxIK+nDV2mPFO4iraDb/Ul7pXAXaQXd7k/aK11zFxEJIY3cRc6gdRtKSe+bE6h2wawZuoQjLaaRu8gZ0i07wlm9swLVHt27i+Li4jbuSMJMI3eRM0RfjpIzSeEu0k5piqW0hsJdpB3S+jbSWgp3kXYo3ks4GuVLYwp3kQ5Oo3xpisJdpINry1E+aKTfUSncRTqReEb5AEfe38ratWsDT8vUG0H7oXAX6UTiGeUDHIr+T+A3A70RtC+Bwt3MCoBHgVTgaXef2+i41R8fD3wG3Oruf0pwryJyhsXzZqA3gvYlZribWSrwBDAWqAQ2mdlKd9/eoOw6YED9z1XAv9f/r4h0Em39RvCjOY8Feu6h/XsEqmuJjvQmE2TkPgSocPedAGa2FJgANAz3CcB/uLsDG8zsfDO70N0/THjHItLhteSNIKignydAfG8Ea9eujeuvjWQLEu59gPcbbFdy6qi8qZo+wGnDvby8PK5P7MMqGo1C94uS3YZIuxXv5wTx2B1H7VfSc+L6MDrZgoS7NbHPW1CDmRUBX/xNc2Tt2rVbA7x++B082PPdedd/nOw22omegM5FHZ2Lk3QuTrosSFGQcK8EGg4tM4E9LajB3RcBiwDMrNTd84M0GXY6FyfpXJykc3GSzsVJZlYapC7Ikr+bgAFmlmVmZwE3ASsb1awEbrY6Q4GDut4uIpI8MUfu7l5rZncAr1A3FXKxu28zs9vqjy8EVlE3DbKCuqmQU9uuZRERiSXQPHd3X0VdgDfct7DB7w7cHudrL4qzPsx0Lk7SuThJ5+IknYuTAp0Lq8tlEREJE91mT0QkhJIa7mb2bTPbZmbHzaxTfhJuZgVmVm5mFWbWae+tZmaLzWyvmXX66bFmdpGZvWFmO+r/fdyZ7J6SxczSzWyjmf1v/bn4ebJ7SiYzSzWzP5vZS7Fqkz1y3wr8K7AuyX0kRYOlHa4DsoHJZpad3K6SZgnQNt9U6XhqgR+7++XAUOD2Tvz/iyPAKHf/GpAHFNTPyOus7gR2BClMari7+w53L09mD0l2YmkHdz8KfLG0Q6fj7uuA/cnuoz1w9w+/WHjP3Q9R94+5T3K7Sg6vc7h+M63+p1N+UGhmmcA3gaeD1Cd75N7ZnW7ZBhEAzKwf8M/AH5PcStLUX4qIAnuB19y9s56L+cD/B44HKW7zcDez1Wa2tYmfTjlCbSTQsg3SOZnZucBy4C53/0ey+0kWdz/m7nnUffN9iJkNSnJLZ5yZXQ/sdffNQR/T5jfrcPcxbf0aHVigZRuk8zGzNOqC/Tl3X5HsftoDdz9gZmuo+2yms33wfjVwo5mNB9KB7mb2W3f/zukeoMsyyRVkaQfpZOpvfvNrYIe7/zLZ/SSTmfUys/Prfz8HGAP8JalNJYG73+vume7ej7qceL25YIfkT4WcaGaVwDDgZTN7JZn9nGnuXgt8sbTDDuA/3X1bcrtKDjN7HngbuMzMKs3s+8nuKYmuBr4LjDKzaP3P+GQ3lSQXAm+YWRl1g6HX3D3mNEDRN1RFREJJl2VEREJI4S4iEkIKdxGREFK4i4iEkMJdRCSEFO4iIiGkcBcRCSGFu4hICP0f+sHibaW8tmcAAAAASUVORK5CYII=\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/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAbd0lEQVR4nO3de3CU9b3H8fc3IRiJ8VII6hg16GA1QsyhIYI4sMqlkYMyTtuRS6ulBzJatXpO6yk6FlvLVBydVvFykLGIpzXIGRXKKCpaRbwhBLuGi2IpYk3RgiAImojB7/ljF7qGJPsk2WU3eT6vmQz7PL/fs/nub0I+eW6/x9wdEREJr5xMFyAiIpmlIBARCTkFgYhIyCkIRERCTkEgIhJyPTJdQEv69OnjJSUlmS5DRKTLWLNmzcfuXtSRbbMyCEpKSqitrc10GSIiXYaZvd/RbXVoSEQk5BQEIiIhpyAQEQm5rDxHICLdw5dffkl9fT2NjY2ZLqXbyM/Pp7i4mLy8vJS9p4JARNKmvr6ewsJCSkpKMLNMl9PluTs7duygvr6efv36pex9dWhIRNKmsbGR3r17KwRSxMzo3bt3yvewFAQiklYKgdRKx3gqCEREQk7nCEJo7ty51NTUBO4/adIkqqur01iRSPcRjUbZunUrY8eObdd2kUiEO++8k4qKijRV1jrtEYRQTU0N0Wg0UN9oNNqu0BAJu2g0ytKlSzNdRrtojyCkGo8+mS1Dbkjeb/P0w1CNSPr88Y9/ZPbs2ezbt49zzz2XH/3oR0ybNo1Vq1axf/9+KisrWbhwIR9//DEzZsygd+/ebNy4keHDh3P//feTk5PDsmXLuOWWW/jiiy84/fTTeeihhzjqqKNYvXo11113HZ999hlHHHEEzz33HDNmzKChoYFXXnmFG2+8kXHjxnHttdeydu1ampqa+OUvf8n48eNpaGhgypQpbNiwgbPOOouGhoaMjZGCQEQOi+uvvz7wnmhQ5eXl3HXXXa22v/322yxcuJBXX32VvLw8fvzjH7Nx40YuueQSbr75ZhoaGvj+97/PgAEDWL58OatWrWLDhg2ceuqpVFVV8cQTTxCJRJg5cybPP/88BQUF3H777fz2t79l+vTpXHbZZSxcuJDBgwfz6aef0qtXL2699VZqa2u59957Abjpppu48MILmTdvHrt27aKyspJRo0bxwAMP0KtXL+rq6qirq2PQoEEpHZv2UBCISLf15z//mTVr1jB48GAAGhoa6Nu3LzNmzGDw4MHk5+cze/bsg/0rKys57bTTAJg4cSKvvPIK+fn5bNiwgWHDhgGwb98+hg4dysaNGznxxBMPvvfRRx/dYg3Lli1jyZIl3HnnnUDsktq///3vrFixgp/85CcAlJWVUVZWlp5BCEBBICKHRVt/uaeLu3PFFVdw2223fW39Rx99xN69e/nyyy9pbGykoKAAOPTSTDPD3Rk9ejQLFiz4WltdXV2gSzndnccff5xvfvObh7Rly6W1OlksIt3WyJEjeeyxx9i2bRsAO3fu5P3336e6uppf//rXTJ48mZ///OcH+69atYr33nuPr776ioULF3L++eczZMgQXn31VTZt2gTA559/zrvvvsuZZ57J1q1bWb16NQB79uyhqamJwsJC9uzZc/A9v/3tb3PPPffg7gD85S9/AWD48OE88sgjAKxbt466urr0D0grtEcgIt1WaWkpM2fOZMyYMXz11Vfk5eUxfvx4evTowaRJk9i/fz/nnXceL7zwAjk5OQwdOpTp06ezdu1ahg8fzqWXXkpOTg7z589n4sSJfPHFFwDMnDmTM844g4ULF3LttdfS0NDAkUceyfPPP88FF1zArFmzKC8v58Ybb+QXv/gF119/PWVlZbg7JSUlPPnkk1x11VVMmTKFsrIyysvLqayszNg4JQ0CM5sHjAO2ufuAFtpvACYnvN9ZQJG77zSzLcAeYD/Q5O6H/wJZEQm1yy67jMsuu6zFttzcXN544w0Ali9fTq9evVi4cOEh/S688MKDf/knGjx4MCtXrjxkffO+DzzwwCF9jjzySB599NFAnyHdghwamg9Utdbo7ne4e7m7lwM3Ai+5+86ELhfE2xUCIiJZKOkegbuvMLOSgO83EViQtJeISJaJRCJEIpFMl5ERKTtZbGa9iO05PJ6w2oFlZrbGzNqco8DMqs2s1sxqt2/fnqqyREQkiVReNXQx8Gqzw0LD3H0QcBFwtZkNb21jd5/r7hXuXlFUVJTCskREpC2pDIIJNDss5O5b4/9uAxYBmTstLiIiLUpJEJjZMcAI4E8J6wrMrPDAa2AMsC4V309ERFInyOWjC4AI0MfM6oFbgDwAd58T73YpsMzdP0vY9HhgUfzOuR5Ajbs/k7rSRaSrKZn+VErfb8usf0/eZ8sWxo0bx7p1X/87dMaMGQwfPpxRo0a1uN3ixYs544wzKC0tTUmt2SzIVUMTA/SZT+wy08R1m4FzOlqYiEg63XrrrW22L168mHHjxqUkCJqamujRI3vv39UUEyLS7e3fv59p06Zx9tlnM2bMGBoaGvjhD3/IY489BsD06dMpLS2lrKyMn/3sZ7z22mssWbKEG264gfLycv72t78RjUYZMmQIZWVlXHrppXzyySdA7OaxsrIyhg4dyg033MCAAbH7bufPn8/3vvc9Lr74YsaMGcPevXsZOXIkgwYNYuDAgfzpT7Ej6Vu2bOHMM89k6tSpDBgwgMmTJ/P8888zbNgw+vfvz6pVq9I+PgoCEen2/vrXv3L11Vezfv16jj32WB5//F9Xue/cuZNFixaxfv166urquPnmmznvvPO45JJLuOOOO4hGo5x++ulcfvnl3H777dTV1TFw4EB+9atfATBlyhTmzJnD66+/Tm5u7te+7+uvv87DDz/MCy+8QH5+PosWLeLNN9/kxRdf5Kc//enB+Yc2bdrEddddR11dHe+88w41NTW88sor3HnnnfzmN79J+/goCESk2+vXrx/l5eUAfOtb32LLli0H244++mjy8/OZOnUqTzzxBL169Tpk+927d7Nr1y5GjBgBwBVXXMGKFSvYtWsXe/bs4bzzzgNij3VNNHr0aL7xjW8AsVlIb7rpJsrKyhg1ahT/+Mc/+Oc//3mwvoEDB5KTk8PZZ5/NyJEjMTMGDhz4tVrTRUEgIt3eEUcccfB1bm4uTU1NB5d79OjBqlWr+M53vsPixYupqmp1Rp1DHPiLvjUHprcGeOSRR9i+fTtr1qwhGo1y/PHH09jYeEh9OTk5B5dzcnK+Vmu6KAhEJNT27t3L7t27GTt2LHfdddfBp6glTid9zDHHcNxxx/Hyyy8D8Ic//IERI0Zw3HHHUVhYeHDiubYmkdu9ezd9+/YlLy+PF198kffffz+9H6wdsvc0toh0O0Eu9zzc9uzZw/jx42lsbMTd+d3vfgfAhAkTmDZtGrNnz+axxx7j4Ycf5sorr+Tzzz/ntNNO46GHHgLg97//PdOmTaOgoIBIJMIxxxzT4veZPHkyF198MRUVFZSXl3PmmWcets+YjCXbtcmEiooKr62tzXQZ3VYkEmHl5h2cMGlW0r4f1UxnyGm9Wb58efoLk27n7bff5qyzzsp0GWm1d+9ejjrqKABmzZrFhx9+yN13353W79nSuJrZmo7O8qw9AhGRTnjqqae47bbbaGpq4tRTT2X+/PmZLqndFAQiIp3Q1oNvugqdLBaRtMrGw89dWTrGU0EgImmTn5/Pjh07FAYp4u7s2LGD/Pz8lL6vDg2JSNoUFxdTX1+PHjaVOvn5+RQXF6f0PRUEIpI2eXl59OvXL9NlSBI6NCQiEnIKAhGRkFMQiIiEnIJARCTkFAQiIiGnIBARCbmkQWBm88xsm5mta6U9Yma7zSwa/5qR0FZlZhvNbJOZTU9l4SIikhpB9gjmA8me1PCyu5fHv24FMLNc4D7gIqAUmGhmnX8KtIiIpFTSIHD3FcDODrx3JbDJ3Te7+z7gUWB8B95HRETSKFXnCIaa2Vtm9rSZnR1fdxLwQUKf+vi6FplZtZnVmlmtbkcXETl8UhEEbwKnuvs5wD3A4vh6a6FvqzNPuftcd69w94qioqIUlCUiIkF0Ogjc/VN33xt/vRTIM7M+xPYATk7oWgxs7ez3ExGR1Op0EJjZCWZm8deV8ffcAawG+ptZPzPrCUwAlnT2+4mISGolnX3UzBYAEaCPmdUDtwB5AO4+B/gucJWZNQENwASPTT7eZGbXAM8CucA8d1+flk8hIiIdljQI3H1ikvZ7gXtbaVsKLO1YaSIicjjozmIRkZBTEIiIhJyCQEQk5BQEIiIhpyAQEQk5BYGISMgpCEREQk5BICIScgoCEZGQUxCIiIScgkBEJOQUBCIiIZd00jnpOkqmPxWo30ebd6S5EhHpShQE3cTcuXP5qOaeQH33bXuPnn37pbkiEekqdGiom6ipqWHftvcC9e3Ztx8FpZH0FiQiXYb2CLqRnn37ccKkWZkuQ0S6GO0RiIiEnIJARCTkFAQiIiGXNAjMbJ6ZbTOzda20TzazuvjXa2Z2TkLbFjNba2ZRM6tNZeEiIpIaQfYI5gNVbbS/B4xw9zLg18DcZu0XuHu5u1d0rEQREUmnpFcNufsKMytpo/21hMWVQHEK6hIRkcMk1ecI/gN4OmHZgWVmtsbMqtva0MyqzazWzGq3b9+e4rJERKQ1KbuPwMwuIBYE5yesHubuW82sL/Ccmb3j7ita2t7d5xI/rFRRUeGpqktERNqWkj0CMysDHgTGu/vBiWzcfWv8323AIqAyFd9PRERSp9NBYGanAE8AP3D3dxPWF5hZ4YHXwBigxSuPREQkc5IeGjKzBUAE6GNm9cAtQB6Au88BZgC9gfvNDKApfoXQ8cCi+LoeQI27P5OGzyAiIp0Q5KqhiUnapwJTW1i/GTjn0C1ERCSb6M5iEZGQ0+yjklQ0GiUSiQTqO2nSJKqr27xSWESyjIJA2lRQGuGzDctZGeCpZgeeh6AgEOlaFATSpsLyKgrL25ph5F8+qpme5mpEJB10jkBEJOQUBCIiIacgEBEJOQWBiEjIKQhEREJOQSAiEnIKAhGRkFMQiIiEnIJARCTkFAQiIiGnIBARCTkFgYhIyCkIRERCTkEgIhJyCgIRkZBLGgRmNs/MtpnZulbazcxmm9kmM6szs0EJbVVmtjHepsnqRUSyUJA9gvlAW08muQjoH/+qBv4HwMxygfvi7aXARDMr7UyxIiKSekmDwN1XADvb6DIe+F+PWQkca2YnApXAJnff7O77gEfjfUVEJIuk4hzBScAHCcv18XWtrW+RmVWbWa2Z1W7fvj0FZYmISBCpCAJrYZ23sb5F7j7X3SvcvaKoqCgFZYmISBCpeHh9PXBywnIxsBXo2cp6ERHJIqnYI1gCXB6/emgIsNvdPwRWA/3NrJ+Z9QQmxPuKiEgWSbpHYGYLgAjQx8zqgVuAPAB3nwMsBcYCm4DPgSnxtiYzuwZ4FsgF5rn7+jR8BhER6YSkQeDuE5O0O3B1K21LiQWFiIhkKd1ZLCIScgoCEZGQUxCIiIScgkBEJOQUBCIiIacgEBEJOQWBiEjIKQhEREJOQSAiEnIKAhGRkFMQiIiEnIJARCTkFAQiIiGnIBARCTkFgYhIyCkIRERCTkEgIhJyqXh4vchB0WiUSCQSqO+kSZOorq5Ob0EikpSCQFKmoDTCwMa1gfpGo1EABYFIFggUBGZWBdxN7CH0D7r7rGbtNwCTE97zLKDI3Xea2RZgD7AfaHL3ihTVLlmmsLyKLVQF6tu4eXqaqxGRoJIGgZnlAvcBo4F6YLWZLXH3DQf6uPsdwB3x/hcD/+nuOxPe5gJ3/zillYuISEoEOVlcCWxy983uvg94FBjfRv+JwIJUFCciIukXJAhOAj5IWK6PrzuEmfUCqoDHE1Y7sMzM1phZqweEzazazGrNrHb79u0ByhIRkVQIEgTWwjpvpe/FwKvNDgsNc/dBwEXA1WY2vKUN3X2uu1e4e0VRUVGAskREJBWCBEE9cHLCcjGwtZW+E2h2WMjdt8b/3QYsInaoSUREskSQIFgN9DezfmbWk9gv+yXNO5nZMcAI4E8J6wrMrPDAa2AMsC4VhYuISGokvWrI3ZvM7BrgWWKXj85z9/VmdmW8fU6866XAMnf/LGHz44FFZnbge9W4+zOp/AAiItI5ge4jcPelwNJm6+Y0W54PzG+2bjNwTqcqFBGRtNJcQyIiIacgEBEJOQWBiEjIKQhEREJOQSAiEnIKAhGRkFMQiIiEnIJARCTkFAQiIiGnIBARCTkFgYhIyCkIRERCTkEgIhJyCgIRkZBTEIiIhJyCQEQk5AI9mEYkHaLRKJFIJFDfSZMmUV1dnd6CREJKQSAZUVAaYWDj2kB9o9EogIJAJE0UBJIRheVVbKEqUN/GzdPTXI1IuAU6R2BmVWa20cw2mdkh/yvNLGJmu80sGv+aEXRbERHJrKR7BGaWC9wHjAbqgdVmtsTdNzTr+rK7j+vgtiIikiFB9ggqgU3uvtnd9wGPAuMDvn9nthURkcMgSBCcBHyQsFwfX9fcUDN7y8yeNrOz27mtiIhkSJCTxdbCOm+2/CZwqrvvNbOxwGKgf8BtY9/ErBqoBjjllFMClCUiIqkQZI+gHjg5YbkY2JrYwd0/dfe98ddLgTwz6xNk24T3mOvuFe5eUVRU1I6PICIinREkCFYD/c2sn5n1BCYASxI7mNkJZmbx15Xx990RZFsREcmspIeG3L3JzK4BngVygXnuvt7Mroy3zwG+C1xlZk1AAzDB3R1ocds0fRYREemAQDeUxQ/3LG22bk7C63uBe4NuKyIi2UOTzomIhJyCQEQk5BQEIiIhpyAQEQk5BYGISMgpCEREQk7PI5AuoT1PMwM90UykPRQEkvUKSiN8tmE5KzfvCNR/37b3AD3RTCQoBYFkvcLyKgrLgz3NDOCjGj3/SKQ9dI5ARCTkFAQiIiGnIBARCTkFgYhIyCkIRERCTkEgIhJyCgIRkZBTEIiIhJxuKJNuqT1TUmg6Cgk7BYF0O+2ZkkLTUYgEDAIzqwLuJvYA+gfdfVaz9snAz+OLe4Gr3P2teNsWYA+wH2hy94rUlC7SsvZMSaHpKEQCBIGZ5QL3AaOBemC1mS1x9w0J3d4DRrj7J2Z2ETAXODeh/QJ3/ziFdYuISIoEOVlcCWxy983uvg94FBif2MHdX3P3T+KLK4Hi1JYpIiLpEiQITgI+SFiuj69rzX8ATycsO7DMzNaYmQ7EiohkmSDnCKyFdd5iR7MLiAXB+Qmrh7n7VjPrCzxnZu+4+4oWtq0GqgFOOeWUAGWJiEgqBNkjqAdOTlguBrY272RmZcCDwHh3P3i5hrtvjf+7DVhE7FDTIdx9rrtXuHtFUVFR8E8gIiKdEiQIVgP9zayfmfUEJgBLEjuY2SnAE8AP3P3dhPUFZlZ44DUwBliXquJFRKTzkh4acvcmM7sGeJbY5aPz3H29mV0Zb58DzAB6A/ebGfzrMtHjgUXxdT2AGnd/Ji2fRKSDdPOZhF2g+wjcfSmwtNm6OQmvpwJTW9huM3BOJ2sUSRvdfCaiO4sl5HTzmYgmnRMRCT0FgYhIyOnQkEg76MSydEcKApGAdGJZuisFgUhAOrEs3ZXOEYiIhJz2CETSROcTpKtQEIikQUFphIGNawP1jUajgM4nSOYoCETSoLC8ii0EO5/QuFnnEySzFAQiWUCHkSSTFAQiGabDSJJpCgKRDNNhJMk0BYFIF9Oew0igQ0mSnIJApAtpz93NAF98sI6XXnqJmpqaQP0VGuGkIBDpQtpzdzPAnugzmhZDklIQiHRj7Z0WQ1cvhZOCQESA9h120iGn7kVBICJA+/Ye2nPI6UBo/GTmPYHee8hpvQP1AwVMqigIRKTdOhIaQaXrRHh7hSlkAgWBmVUBdwO5wIPuPqtZu8XbxwKfAz909zeDbCsi3Vt7T3AH1Z69kvZKd8hkm6RBYGa5wH3AaKAeWG1mS9x9Q0K3i4D+8a9zgf8Bzg24rYhIu6UrYCC9IZONguwRVAKb3H0zgJk9CowHEn+Zjwf+190dWGlmx5rZiUBJgG0PsXHjxnbdMCPxqQeOPjnTZYh0C+kMmXR5//ZxHd42SBCcBHyQsFxP7K/+ZH1OCrgtAGZWDRw4IPfFSy+9tC5AbZnUB/g400V8ze7dLf0wZF+dLVOdqaU6U6sr1PnNjm4YJAishXUesE+QbWMr3ecCcwHMrNbdKwLUljFdoUZQnammOlNLdaaOmdV2dNsgQVAPJB5zKAa2BuzTM8C2IiKSQUGeWbwa6G9m/cysJzABWNKszxLgcosZAux29w8DbisiIhmUdI/A3ZvM7BrgWWKXgM5z9/VmdmW8fQ6wlNilo5uIXT46pa1tA9Q1tyMf5jDrCjWC6kw11ZlaqjN1OlyjxS70ERGRsApyaEhERLoxBYGISMhlLAjMrMrMNprZJjM75Pl78RPPs+PtdWY2KEvrjJjZbjOLxr9mZKDGeWa2zcxavPcii8YyWZ0ZH8t4HSeb2Ytm9raZrTez61rok/ExDVhnRsfUzPLNbJWZvRWv8Vct9MmGsQxSZ1b8fMZryTWzv5jZky20tX883f2wfxE7cfw34DRil5i+BZQ26zMWeJrYvQhDgDeytM4I8GQmxjGhhuHAIGBdK+0ZH8uAdWZ8LON1nAgMir8uBN7N0p/PIHVmdEzj43NU/HUe8AYwJAvHMkidWfHzGa/lv4CalurpyHhmao/g4LQV7r4PODD1RKKD01a4+0rgwLQV2VZnxrn7CmBnG12yYSyD1JkV3P1Dj0+a6O57gLeJ3SWfKONjGrDOjIqPz974Yl78q/kVKtkwlkHqzApmVgz8O/BgK13aPZ6ZCoLWpqRob590C1rD0Pgu5dNmdvbhKa1dsmEsg8qqsTSzEuDfiP2FmCirxrSNOiHDYxo/jBEFtgHPuXtWjmWAOiE7fj7vAv4b+KqV9naPZ6aCoDPTVhxOQWp4EzjV3c8B7gEWp7uoDsiGsQwiq8bSzI4CHgeud/dPmze3sElGxjRJnRkfU3ff7+7lxGYWqDSzAc26ZMVYBqgz42NpZuOAbe6+pq1uLaxrczwzFQSdmbbicEpag7t/emCX0t2XAnlm1ufwlRhINoxlUtk0lmaWR+yX6yPu/kQLXbJiTJPVmU1j6u67gOVA82k9s2IsD2itziwZy2HAJWa2hdih6gvN7I/N+rR7PDMVBJ2ZtiKr6jSzE8zM4q8riY1ptk1ing1jmVS2jGW8ht8Db7v7b1vplvExDVJnpsfUzIrM7Nj46yOBUcA7zbplw1gmrTPTYwng7je6e7G7lxD7ffSCu3+/Wbd2j2dGHlXpnZi2Igvr/C5wlZk1AQ3ABI+fuj9czGwBsSsa+phZPXALsZNdWTOWAevM+FjGDQN+AKyNHzMGuAk4JaHWbBjTIHVmekxPBB622EOqcoD/c/cns+3/esA6Mz2WrerseGqKCRGRkNOdxSIiIacgEBEJOQWBiEjIKQhEREJOQSAiEnIKAhGRkFMQiIiE3P8DaxEQfa+3TIsAAAAASUVORK5CYII=\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/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAauklEQVR4nO3df5RVdb3/8eeLEZokUNNQF4MxejEaYZhoGAVcOCYSmsp1dV38Ks2uTlRmrm/5jVzGd6V+S1auMst7aVaR3pLgLhWbpZTmDxb+Ihhq7ggoNiHmhAZiIiiIyPv7xxz4HsczzJ5hhn1gvx5rzeLs/fl8znmfD6x5sT97n30UEZiZWXb1SbsAMzNLl4PAzCzjHARmZhnnIDAzyzgHgZlZxh2RdgGFHHfccTF06NC0yzAzO2SsWrXq1Yj4SHfGFmUQDB06lMbGxrTLMDM7ZEh6sbtjvTRkZpZxDgIzs4xzEJiZZVxRniMo5J133qG1tZWdO3emXcpho7S0lLKyMvr27Zt2KWaWokMmCFpbWxkwYABDhw5FUtrlHPIigi1bttDa2kp5eXna5ZhZig6ZpaGdO3dy7LHHOgR6iCSOPfZYH2GZWbIgkDRZ0jpJLZJm76ffGEnvSvq3ro5NWMeBDLd2PJ9mBgmCQFIJcDtwHlABTJdU0UG/ucCDXR1rZmbpSXJEUAO0RMT6iNgFLASmFOj3NeAeYFM3xloBTU1NLFmypMvjamtr/YE8M0ssycniwcBLedutwOn5HSQNBi4GPgWM6cpY61hTUxONjY2cf/75aZdi1uPq6+u5+qafpF2GkeyIoNBCcvuvNbsV+FZEvNuNsW0dpTpJjZIaN2/enKCsg+/Xv/41NTU1VFVV8aUvfYk//vGPVFZWsnPnTt58801OO+00Vq9ezdKlS5kwYQIXX3wxFRUVzJo1iz179gDw0EMPMXbsWEaPHs0ll1zC9u3bAVi5ciXjxo1j1KhR1NTUsHXrVubMmcOiRYuoqqpi0aJFvPnmm3zxi19kzJgxfOITn+C3v/0tADt27GDatGlUVlYydepUduzYkdocmSW1YMECdm16Ie0yjGRHBK3AkLztMmBjuz7VwMLcycfjgPMl7U44FoCIqAfqAaqrq/f7/ZnXXHMNTU1NCUpPrqqqiltvvbXD9meffZZFixbx5JNP0rdvX77yla+wbt06LrroIq6//np27NjB5z73OUaMGMHSpUtZsWIFa9eu5aMf/SiTJ0/m3nvvpba2lptuuomHH36Y/v37M3fuXH74wx8ye/Zspk6dyqJFixgzZgxvvPEGRx55JDfccAONjY389Kc/BeC6667jU5/6FPPnz+f111+npqaGiRMn8rOf/YwjjzyS5uZmmpubGT16dI/OjVlv6TeonBNm3Jx2GYeFF+de0O2xSYJgJTBMUjnwd2AaMCO/Q0TsuxBd0h3A/RFxn6QjOht7qHjkkUdYtWoVY8a0rXzt2LGDQYMGMWfOHMaMGUNpaSm33Xbbvv41NTWcfPLJAEyfPp0nnniC0tJS1q5dy/jx4wHYtWsXY8eOZd26dZx44on7nnvgwIEFa3jooYdoaGjglltuAdouqf3b3/7GsmXLuPrqqwGorKyksrKydybBzA5LnQZBROyWdBVtVwOVAPMjYo2kWbn2eV0de6BF7+9/7r0lIrjsssv4/ve//579r7zyCtu3b+edd95h586d9O/fH3j/pZmSiAjOPfdcfvOb37ynrbm5OdGlnBHBPffcw8c+9rH3tflSUDPrrkSfI4iIJRFxakScEhH/N7dvXqEQiIgvRMTd+xt7KDrnnHO4++672bSp7aKo1157jRdffJG6ujpuvPFGZs6cybe+9a19/VesWMELL7zAnj17WLRoEWeeeSZnnHEGTz75JC0tLQC89dZbPP/88wwfPpyNGzeycuVKALZt28bu3bsZMGAA27Zt2/ecn/70p/nJT35CRNvK2Z///GcAJkyYwF133QXA6tWraW5u7v0JMbPDxiFzi4m0VVRUcNNNNzFp0iT27NlD3759mTJlCkcccQQzZszg3XffZdy4cTz66KP06dOHsWPHMnv2bJ555pl9J4779OnDHXfcwfTp03n77bcBuOmmmzj11FNZtGgRX/va19ixYwcf/OAHefjhhzn77LO5+eabqaqq4tvf/jbf+c53uOaaa6isrCQiGDp0KPfffz9f/vKXufzyy6msrKSqqoqampqUZ8vMDiXa+7/LYlJdXR3tr4N/9tln+fjHP55SRV2zdOlSbrnlFu6///60S+nUoTSvdnipra1l+fotPlncQ16ce8GqiKjuzthD5l5DZmbWO7w01Atqa2upra1Nuwwzs0R8RGBmlnEOAjOzjHMQmJllnIPAzCzjDtmTxUNnP9Cjz7fh5s/sv33DBi644AJWr179nv1z5sxhwoQJTJw4seC4++67j1NPPZWKCn8Ng5kVJx8RHKAbbrihwxCAtiBYu3Ztj7zW7t27e+R5zMzyOQi64N133+XKK6/ktNNOY9KkSezYsYMvfOEL3H132x01Zs+eTUVFBZWVlXzzm9/kqaeeoqGhgWuvvZaqqir++te/0tTUxBlnnEFlZSUXX3wx//znP4G221BXVlYyduxYrr32WkaMGAHAHXfcwSWXXMKFF17IpEmT2L59O+eccw6jR49m5MiR+25FvWHDBoYPH84VV1zBiBEjmDlzJg8//DDjx49n2LBhrFixIp1JM7Oi5yDogr/85S989atfZc2aNRx99NHcc889+9pee+01Fi9ezJo1a2hubub6669n3LhxXHTRRfzgBz+gqamJU045hUsvvZS5c+fS3NzMyJEj+e53vwvA5Zdfzrx583j66acpKSl5z+s+/fTT3HnnnTz66KOUlpayePFi/vSnP/HYY4/xjW98Y9+9h1paWvj6179Oc3Mzzz33HAsWLOCJJ57glltu4Xvf+97BmygzO6Q4CLqgvLycqqoqAD75yU+yYcOGfW0DBw6ktLSUK664gnvvvZcjjzzyfeO3bt3K66+/zllnnQXAZZddxrJly3j99dfZtm0b48aNA2DGjPfeqfvcc8/lwx/+MNB2B9LrrruOyspKJk6cyN///nf+8Y9/7Ktv5MiR9OnTh9NOO41zzjkHSYwcOfI9tZqZ5XMQdMEHPvCBfY9LSkres2Z/xBFHsGLFCj772c9y3333MXny5MTP29n9nvbe2hrgrrvuYvPmzaxatYqmpiaOP/54du7c+b76+vTps2+7T58+Pr9gZh1yEPSQ7du3s3XrVs4//3xuvfXWfd+gln8r6aOOOopjjjmGxx9/HIBf/epXnHXWWRxzzDEMGDCA5cuXA7Bw4cIOX2fr1q0MGjSIvn378thjj/Hiiy/27hszs8PeIXv5aGeXex5s27ZtY8qUKezcuZOI4Ec/+hEA06ZN48orr+S2227j7rvv5s4772TWrFm89dZbnHzyyfzyl78E4Be/+AVXXnkl/fv3p7a2lqOOOqrg68ycOZMLL7yQ6upqqqqqGD58+EF7j2Z2ePJtqIvE9u3b+dCHPgTAzTffzMsvv8yPf/zjXn/dw31erXj5NtQ9q9dvQy1psqR1klokzS7QPkVSs6QmSY2Szsxr2yDpmb1t3SkyCx544AGqqqoYMWIEjz/+ONdff33aJZlZRnS6NCSpBLgdOBdoBVZKaoiI/E9JPQI0RERIqgT+G8hfszg7Il7twboPO1OnTmXq1Klpl2FmGZTkiKAGaImI9RGxC1gITMnvEBHb4/+vMfUHemW9qRiXsQ5lnk8zg2RBMBh4KW+7NbfvPSRdLOk54AHgi3lNATwkaZWkuo5eRFJdblmpcfPmze9rLy0tZcuWLf7l1UMigi1btlBaWpp2KWaWsiRXDanAvvf9No6IxcBiSROAG4G9N+AZHxEbJQ0C/iDpuYhYVmB8PVAPbSeL27eXlZXR2tpKoZCw7iktLaWsrCztMswsZUmCoBUYkrddBmzsqHNELJN0iqTjIuLViNiY279J0mLalpreFwSd6du3L+Xl5V0dZmZmnUiyNLQSGCapXFI/YBrQkN9B0r9IUu7xaKAfsEVSf0kDcvv7A5OA997H2czMUtXpEUFE7JZ0FfAgUALMj4g1kmbl2ucBnwUulfQOsAOYmruC6Hjalov2vtaCiPh9L70XMzPrhkSfLI6IJcCSdvvm5T2eC8wtMG49MOoAazQzs17kew2ZmWWcg8DMLOMcBGZmGecgMDPLOAeBmVnGOQjMzDLOQWBmlnEOAjOzjHMQmJllnIPAzCzjHARmZhnnIDAzyzgHgZlZxjkIzMwyzkFgZpZxDgIzs4xzEJiZZVyiIJA0WdI6SS2SZhdonyKpWVKTpEZJZyYda2Zm6eo0CCSVALcD5wEVwHRJFe26PQKMiogq4IvAz7sw1szMUpTkO4trgJbc9w8jaSEwBVi7t0NEbM/r3x+IpGOt99TX17NgwYK0yziszJgxg7q6urTLMOtRSYJgMPBS3nYrcHr7TpIuBr4PDAI+05WxufF1QB3ASSedlKAs68yCBQtYtryRfoPK0y7lsLBr0wsADgI77CQJAhXYF+/bEbEYWCxpAnAjMDHp2Nz4eqAeoLq6umAf67p+g8o5YcbNaZdxWHhlgU9x2eEpycniVmBI3nYZsLGjzhGxDDhF0nFdHWtmZgdfkiBYCQyTVC6pHzANaMjvIOlfJCn3eDTQD9iSZKyZmaWr06WhiNgt6SrgQaAEmB8RayTNyrXPAz4LXCrpHWAHMDUiAig4tpfei5mZdUOScwRExBJgSbt98/IezwXmJh1rZmbFw58sNjPLOAeBmVnGOQjMzDLOQWBmlnEOAjOzjHMQmJllnIPAzCzjHARmZhnnIDAzyzgHgZlZxjkIzMwyzkFgZpZxDgIzs4xzEJiZZZyDwMws4xwEZmYZ5yAwM8u4REEgabKkdZJaJM0u0D5TUnPu5ylJo/LaNkh6RlKTpMaeLN7MzA5cp19VKakEuB04F2gFVkpqiIi1ed1eAM6KiH9KOg+oB07Paz87Il7twbrNzKyHJDkiqAFaImJ9ROwCFgJT8jtExFMR8c/c5nKgrGfLNDOz3pIkCAYDL+Vtt+b2deTfgd/lbQfwkKRVkuo6GiSpTlKjpMbNmzcnKMvMzHpCp0tDgArsi4IdpbNpC4Iz83aPj4iNkgYBf5D0XEQse98TRtTTtqREdXV1wec3M7Oel+SIoBUYkrddBmxs30lSJfBzYEpEbNm7PyI25v7cBCymbanJzMyKRJIgWAkMk1QuqR8wDWjI7yDpJOBe4PMR8Xze/v6SBux9DEwCVvdU8WZmduA6XRqKiN2SrgIeBEqA+RGxRtKsXPs8YA5wLPAfkgB2R0Q1cDywOLfvCGBBRPy+V96JmZl1S5JzBETEEmBJu33z8h5fAVxRYNx6YFT7/WZmVjz8yWIzs4xzEJiZZZyDwMws4xwEZmYZ5yAwM8s4B4GZWcY5CMzMMs5BYGaWcQ4CM7OMcxCYmWWcg8DMLOMcBGZmGecgMDPLOAeBmVnGOQjMzDLOQWBmlnEOAjOzjEsUBJImS1onqUXS7ALtMyU1536ekjQq6VgzM0tXp0EgqQS4HTgPqACmS6po1+0F4KyIqARuBOq7MNbMzFKU5DuLa4CW3PcPI2khMAVYu7dDRDyV1385UJZ0rNmhpKmpidra2rTLOCw0NTXBwCFpl2EkC4LBwEt5263A6fvp/+/A77o6VlIdUAdw0kknJSjL7ODqX1HLm2uXsnz9lrRLOTwMHEL/itq0qzCSBYEK7IuCHaWzaQuCM7s6NiLqyS0pVVdXF+xjlqYBVZMZUDU57TLMelySIGgF8o/fyoCN7TtJqgR+DpwXEVu6MtbMzNKT5KqhlcAwSeWS+gHTgIb8DpJOAu4FPh8Rz3dlrJmZpavTI4KI2C3pKuBBoASYHxFrJM3Ktc8D5gDHAv8hCWB3RFR3NLaX3ouZmXVDkqUhImIJsKTdvnl5j68Arkg61szMioc/WWxmlnEOAjOzjHMQmJllnIPAzCzjHARmZhnnIDAzyzgHgZlZxjkIzMwyzkFgZpZxDgIzs4xzEJiZZZyDwMws4xwEZmYZ5yAwM8s4B4GZWcY5CMzMMi5REEiaLGmdpBZJswu0D5f0tKS3JX2zXdsGSc9IapLU2FOFm5lZz+j0G8oklQC3A+fS9mX0KyU1RMTavG6vAVcD/9rB05wdEa8eYK1mZtYLkhwR1AAtEbE+InYBC4Ep+R0iYlNErATe6YUazcysFyUJgsHAS3nbrbl9SQXwkKRVkuo66iSpTlKjpMbNmzd34enNzOxAJAkCFdgXXXiN8RExGjgP+KqkCYU6RUR9RFRHRPVHPvKRLjy9mZkdiCRB0AoMydsuAzYmfYGI2Jj7cxOwmLalJjMzKxJJgmAlMExSuaR+wDSgIcmTS+ovacDex8AkYHV3izUzs57X6VVDEbFb0lXAg0AJMD8i1kialWufJ+kEoBEYCOyRdA1QARwHLJa097UWRMTve+WdmJlZt3QaBAARsQRY0m7fvLzHr9C2ZNTeG8CoAynQzMx6lz9ZbGaWcQ4CM7OMcxCYmWWcg8DMLOMcBGZmGecgMDPLOAeBmVnGOQjMzDLOQWBmlnEOAjOzjHMQmJllnIPAzCzjHARmZhnnIDAzyzgHgZlZxjkIzMwyzkFgZpZxiYJA0mRJ6yS1SJpdoH24pKclvS3pm10Za2Zm6eo0CCSVALcD59H2PcTTJVW06/YacDVwSzfGmplZipJ8Z3EN0BIR6wEkLQSmAGv3doiITcAmSZ/p6thC1q1bR21tbdL3YB1oamqCgUPSLsPMilySIBgMvJS33QqcnvD5E4+VVAfUAVDSl+XrtyR8CevQwCH0r6hNuwozK3JJgkAF9kXC5088NiLqgXqAD5w4LE6YcXPClzAzswOR5GRxK5C/vlAGbEz4/Acy1szMDoIkQbASGCapXFI/YBrQkPD5D2SsmZkdBJ0uDUXEbklXAQ8CJcD8iFgjaVaufZ6kE4BGYCCwR9I1QEVEvFFobC+9FzMz64Yk5wiIiCXAknb75uU9foW2ZZ9EY83MrHj4k8VmZhnnIDAzyzgHgZlZxjkIzMwyzkFgZpZxDgIzs4xzEJiZZZyDwMws4xwEZmYZ5yAwM8s4B4GZWcY5CMzMMs5BYGaWcQ4CM7OMcxCYmWWcg8DMLOMcBGZmGZcoCCRNlrROUouk2QXaJem2XHuzpNF5bRskPSOpSVJjTxZvZmYHrtOvqpRUAtwOnAu0AislNUTE2rxu5wHDcj+nA/+Z+3OvsyPi1R6r2szMekySI4IaoCUi1kfELmAhMKVdnynAf0Wb5cDRkk7s4VrNzKwXJAmCwcBLedutuX1J+wTwkKRVkuo6ehFJdZIaJTW++9bWBGWZmVlP6HRpCFCBfdGFPuMjYqOkQcAfJD0XEcve1zmiHqgH+MCJw9o/v5mZ9ZIkRwStwJC87TJgY9I+EbH3z03AYtqWmszMrEgkCYKVwDBJ5ZL6AdOAhnZ9GoBLc1cPnQFsjYiXJfWXNABAUn9gErC6B+s3M7MD1OnSUETslnQV8CBQAsyPiDWSZuXa5wFLgPOBFuAt4PLc8OOBxZL2vtaCiPh9j78LMzPrtiTnCIiIJbT9ss/fNy/vcQBfLTBuPTDqAGs0M7Ne5E8Wm5llnIPAzCzjHARmZhnnIDAzyzgHgZlZxjkIzMwyzkFgZpZxDgIzs4xzEJiZZZyDwMws4xwEZmYZ5yAwM8s4B4GZWcY5CMzMMs5BYGaWcQ4CM7OMcxCYmWVcoiCQNFnSOkktkmYXaJek23LtzZJGJx1rZmbp6jQIJJUAtwPnARXAdEkV7bqdBwzL/dQB/9mFsWZmlqIkRwQ1QEtErI+IXcBCYEq7PlOA/4o2y4GjJZ2YcKyZmaUoyZfXDwZeyttuBU5P0GdwwrEASKqj7WgC4O0X516wOkFtaToOeDXtIhJwnT3LdfYs19lzPtbdgUmCQAX2RcI+Sca27YyoB+oBJDVGRHWC2lJzKNQIrrOnuc6e5Tp7jqTG7o5NEgStwJC87TJgY8I+/RKMNTOzFCU5R7ASGCapXFI/YBrQ0K5PA3Bp7uqhM4CtEfFywrFmZpaiTo8IImK3pKuAB4ESYH5ErJE0K9c+D1gCnA+0AG8Bl+9vbIK66rvzZg6yQ6FGcJ09zXX2LNfZc7pdoyIKLtmbmVlG+JPFZmYZ5yAwM8u41IIgwW0raiVtldSU+5mTUp3zJW2SVPBzDfu7vUYR1VgsczlE0mOSnpW0RtLXC/QphvlMUmfqcyqpVNIKSf+Tq/O7BfqkOp8Ja0x9LvNqKZH0Z0n3F2hL/d9mXi37q7Pr8xkRB/2HthPHfwVOpu0S0/8BKtr1qQXuT6O+dnVMAEYDqztoPx/4HW2fmTgD+GMR1lgsc3kiMDr3eADwfIG/92KYzyR1pj6nuTn6UO5xX+CPwBnFNJ8Ja0x9LvNq+V/AgkL1pD2XXaizy/OZ1hHBIXPriYhYBry2ny4d3V7joElQY1GIiJcj4k+5x9uAZ2n79Hm+YpjPJHWmLjdH23ObfXM/7a/+SHU+E9ZYFCSVAZ8Bft5Bl9T/bUKiOrssrSDo6JYU7Y3NHVL+TtJpB6e0Lkv6XtJWVHMpaSjwCdr+h5ivqOZzP3VCEcxpbomgCdgE/CEiim4+E9QIRTCXwK3A/wb2dNCe+lzm3Mr+64QuzmdaQZDk1hN/Aj4aEaOAnwD39XZR3ZT4NhopKqq5lPQh4B7gmoh4o31zgSGpzGcndRbFnEbEuxFRRdun9mskjWjXJfX5TFBj6nMp6QJgU0Ss2l+3AvsO6lwmrLPL85lWEHR624qIeGPvIWVELAH6Sjru4JWYWJJbcKSqmOZSUl/afrneFRH3FuhSFPPZWZ3FNKe5Gl4HlgKT2zUVxXxCxzUWyVyOBy6StIG2pepPSfp1uz7FMJed1tmd+UwrCDq99YSkEyQp97iGtlq3HPRKO9fR7TWKRrHMZa6GXwDPRsQPO+iW+nwmqbMY5lTSRyQdnXv8QWAi8Fy7bqnOZ5Iai2EuI+LbEVEWEUNp+330aER8rl231P9tJqmzO/OZ5KZzPS6S3bbi34AvS9oN7ACmRe6U+MEk6Te0nYU/TlIr8H9oO+G1t86Ct9coshqLYi5p+9/M54FncmvGANcBJ+XVmvp8kqzOYpjTE4E71fYFUH2A/46I+5Xg9i9FVmMxzGVBRTaXHTrQ+fQtJszMMs6fLDYzyzgHgZlZxjkIzMwyzkFgZpZxDgIzs4xzEJiZZZyDwMws4/4f1QSTbeiJDgsAAAAASUVORK5CYII=\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/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAX8UlEQVR4nO3dfXBV9b3v8fc3IRBBQAulMoQj0YvFCGkOhQjogbQ8NHLQHKbt8NSjpVcYbLE60zqNjqWtciqMTEvx4VLmFPG2UjiDSDNIr5Yql+LDIcGm4cnYVEUjWhQKggQx+D1/ZMts445ZSfZmZ+f3ec1kZj389trf3/4xH1Z+e60Vc3dERKTry0p3ASIicm4o8EVEAqHAFxEJhAJfRCQQCnwRkUB0S9cb9+/f34cMGZKutxcRyUi7du16x90/257Xpi3whwwZQlVVVbreXkQkI5nZgfa+VlM6IiKBUOCLiARCgS8iEoi0zeGLSNfwwQcfUF9fz6lTp9JdSpeSm5tLXl4eOTk5STumAl9EOqS+vp7evXszZMgQzCzd5XQJ7s7hw4epr68nPz8/acfVlI6IdMipU6fo16+fwj6JzIx+/fol/bemSIFvZqVmVmtmdWZWnmD/bWZWHfvZY2ZnzOwzSa1URDothX3ypeIzbTXwzSwbeAC4BigAZplZQXwbd7/X3YvcvQi4Hfj/7n4k6dWKiEi7RZnDLwbq3P1lADNbB5QB+1poPwv4bWsHra2tpaSkJGKZ0tnMnj2b+fPnp7sMkbSprq7m4MGDTJ06tU2vKykpYdmyZYwaNSpFlbUsSuAPAl6PW68HrkzU0Mx6AqXAwhb2zweaUiI7h+dfPtyWWqWTOH3oFQAFvgSturqaqqqqNgd+OkUJ/EQTSS39maxrgWdams5x91XAKoAeA4f6RbOXRCpSOpe31n7iaxyRtPvNb37DihUrOH36NFdeeSXf+ta3mDdvHjt37uTMmTMUFxezfv163nnnHRYtWkS/fv2ora1l/PjxPPjgg2RlZfHkk0/yox/9iPfff59LL72Uhx56iPPPP5/KykpuueUW3nvvPXr06MEf/vAHFi1aRENDAzt27OD2229n2rRp3HzzzezevZvGxkZ+/OMfU1ZWRkNDA3PnzmXfvn1cfvnlNDQ0pO0zihL49cDguPU84GALbWcSYTpHRLqmW2+9lerq6qQes6ioiOXLl39qm/3797N+/XqeeeYZcnJy+Pa3v01tbS3XXXcdd955Jw0NDXzjG99g+PDhbNu2jZ07d7Jv3z4uvvhiSktL2bhxIyUlJSxevJitW7fSq1cvli5dys9+9jPKy8uZMWMG69evZ/To0bz77rv07NmTu+66i6qqKu6//34A7rjjDr785S+zevVqjh49SnFxMZMmTeKXv/wlPXv2pKamhpqaGkaOHJnUz6ctogR+JTDUzPKBN2gK9dnNG5lZX2AC8I2kVigi0oo//vGP7Nq1i9GjRwPQ0NDAgAEDWLRoEaNHjyY3N5cVK1acbV9cXMwll1wCwKxZs9ixYwe5ubns27ePq666CoDTp08zduxYamtrGThw4Nlj9+nTJ2ENTz75JBUVFSxbtgxoulz1tddeY/v27Xz3u98FoLCwkMLCwtR8CBG0Gvju3mhmC4EngGxgtbvvNbMFsf0rY02nA0+6+3spq1ZEOrXWzsRTxd254YYbuOeeez62/a233uLEiRN88MEHnDp1il69egGfvOTRzHB3Jk+ezG9/+/FJipqamkiXSLo7jz76KJ///Oc/sa+zXLYa6Tp8d9/i7pe5+6Xu/h+xbSvjwh53X+PuM1NVqIhISyZOnMiGDRs4dOgQAEeOHOHAgQPMnz+fu+++mzlz5vCDH/zgbPudO3fyyiuv8OGHH7J+/XquvvpqxowZwzPPPENdXR0AJ0+e5KWXXmLYsGEcPHiQyspKAI4fP05jYyO9e/fm+PHjZ4/5la98hfvuuw/3pq84//znPwMwfvx4HnnkEQD27NlDTU1N6j+QFujRCiKS8QoKCli8eDFTpkzhww8/JCcnh7KyMrp168bs2bM5c+YM48aN46mnniIrK4uxY8dSXl7O7t27GT9+PNOnTycrK4s1a9Ywa9Ys3n//fQAWL17MZZddxvr167n55ptpaGjgvPPOY+vWrXzpS19iyZIlFBUVcfvtt/PDH/6QW2+9lcLCQtydIUOGsHnzZm666Sbmzp1LYWEhRUVFFBcXp+1zso/+NzrXegwc6gNvWJ6W95aOeWttOWMu6ce2bdvSXYp0Avv37+fyyy9PdxmRbdu2jWXLlrF58+Z0l9KqRJ+tme1y93ZdxK9n6YiIBEJTOiISlJKSkmDv8tcZvohIIBT4IiKBUOCLiARCgS8iEgh9aSsiSTWk/PGkHu/VJf/aeptXX2XatGns2bPnY9sXLVrE+PHjmTRpUsLXbdq0icsuu4yCgoKE+7saneGLSJd11113tRj20BT4+/a19Kc92qaxsTEpx0klBb6IdAlnzpxh3rx5XHHFFUyZMoWGhga++c1vsmHDBgDKy8spKCigsLCQ73//+zz77LNUVFRw2223UVRUxN/+9jeqq6sZM2YMhYWFTJ8+nX/84x8AVFZWUlhYyNixY7ntttsYPnw4AGvWrOHrX/861157LVOmTOHEiRNMnDiRkSNHMmLECH73u98BTb+BDBs2jBtvvJHhw4czZ84ctm7dylVXXcXQoUPZuXPnOfmMFPgi0iX89a9/5Tvf+Q579+7lggsu4NFHHz2778iRIzz22GPs3buXmpoa7rzzTsaNG8d1113HvffeS3V1NZdeeinXX389S5cupaamhhEjRvCTn/wEgLlz57Jy5Uqee+45srOzP/a+zz33HA8//DBPPfUUubm5PPbYY7zwwgs8/fTTfO973zv7bJ26ujpuueUWampqePHFF1m7di07duxg2bJl/PSnPz0nn5ECX0S6hPz8fIqKigD44he/yKuvvnp2X58+fcjNzeXGG29k48aN9OzZ8xOvP3bsGEePHmXChAkA3HDDDWzfvp2jR49y/Phxxo0bBzT9ec94kydP5jOf+QzQ9MTMO+64g8LCQiZNmsQbb7zB3//+97P1jRgxgqysLK644gomTpyImTFixIiP1ZpKCnwR6RJ69Ohxdjk7O/tjc+rdunVj586dfPWrX2XTpk2UlpZGPm5rzxv76JHLAI888ghvv/02u3btorq6ms997nOcOnXqE/VlZWWdXc/Kyjpn8/8KfBHp8k6cOMGxY8eYOnUqy5cvP/tXueIfcdy3b18uvPBC/vSnPwHw61//mgkTJnDhhRfSu3dvnn/+eQDWrVvX4vscO3aMAQMGkJOTw9NPP82BAwdS27E20mWZIpJUUS6jPNeOHz9OWVkZp06dwt35+c9/DsDMmTOZN28eK1asYMOGDTz88MMsWLCAkydPcskll/DQQw8B8Ktf/Yp58+bRq1cvSkpK6Nu3b8L3mTNnDtdeey2jRo2iqKiIYcOGnbM+RqHHI0ub6fHIEi/THo/cHidOnOD8888HYMmSJbz55pv84he/SPn7JvvxyDrDFxFpxeOPP84999xDY2MjF198MWvWrEl3Se2iwBcRacWMGTOYMWNGusvoMH1pKyIdlq6p4a4sFZ+pAl9EOiQ3N5fDhw8r9JPI3Tl8+DC5ublJPW6kKR0zKwV+AWQD/+nuSxK0KQGWAznAO+4+IWlVikinlZeXR319PW+//Xa6S+lScnNzycvLS+oxWw18M8sGHgAmA/VApZlVuPu+uDYXAA8Cpe7+mpkNSGqVItJp5eTkkJ+fn+4yJIIoUzrFQJ27v+zup4F1QFmzNrOBje7+GoC7H0pumSIi0lFRAn8Q8Hrcen1sW7zLgAvNbJuZ7TKz6xMdyMzmm1mVmVWdOXmsfRWLiEi7RJnDtwTbmn870w34IjAROA94zsyed/eXPvYi91XAKmi68art5YqISHtFCfx6YHDceh5wMEGbd9z9PeA9M9sOfAF4CRER6RSiTOlUAkPNLN/MugMzgYpmbX4H/IuZdTOznsCVwP7klioiIh3R6hm+uzea2ULgCZouy1zt7nvNbEFs/0p3329m/w+oAT6k6dLNPS0fVUREzrVI1+G7+xZgS7NtK5ut3wvcm7zSREQkmXSnrYhIIBT4IiKBUOCLiARCgS8iEggFvohIIBT4IiKBUOCLiARCgS8iEggFvohIIBT4IiKBUOCLiARCgS8iEggFvohIIBT4IiKBUOCLiARCgS8iEggFvohIIBT4IiKBUOCLiARCgS8iEohIgW9mpWZWa2Z1ZlaeYH+JmR0zs+rYz6LklyoiIh3RrbUGZpYNPABMBuqBSjOrcPd9zZr+yd2npaBGERFJgihn+MVAnbu/7O6ngXVAWWrLEhGRZIsS+IOA1+PW62PbmhtrZn8xs9+b2RWJDmRm882sysyqzpw81o5yRUSkvaIEviXY5s3WXwAudvcvAPcBmxIdyN1Xufsodx+V3bNvmwoVEZGOiRL49cDguPU84GB8A3d/191PxJa3ADlm1j9pVYqISIdFCfxKYKiZ5ZtZd2AmUBHfwMwuMjOLLRfHjns42cWKiEj7tXqVjrs3mtlC4AkgG1jt7nvNbEFs/0rga8BNZtYINAAz3b35tI+IiKRRq4EPZ6dptjTbtjJu+X7g/uSWJiIiyaQ7bUVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAhEp8M2s1MxqzazOzMo/pd1oMztjZl9LXokiIpIMrQa+mWUDDwDXAAXALDMraKHdUuCJZBcpIiIdF+UMvxioc/eX3f00sA4oS9DuZuBR4FAS6xMRkSSJEviDgNfj1utj284ys0HAdGDlpx3IzOabWZWZVZ05eayttYqISAdECXxLsM2brS8HfuDuZz7tQO6+yt1Hufuo7J59I5YoIiLJ0C1Cm3pgcNx6HnCwWZtRwDozA+gPTDWzRnfflIwiRUSk46IEfiUw1MzygTeAmcDs+Abunv/RspmtATYr7EVEOpdWA9/dG81sIU1X32QDq919r5ktiO3/1Hl7ERHpHKKc4ePuW4AtzbYlDHp3/2bHyxIRkWTTnbYiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBEKBLyISCAW+iEggFPgiIoFQ4IuIBEKBLyISCAW+iEggIgW+mZWaWa2Z1ZlZeYL9ZWZWY2bVZlZlZlcnv1QREemIbq01MLNs4AFgMlAPVJpZhbvvi2v2R6DC3d3MCoH/AoalomAREWmfKGf4xUCdu7/s7qeBdUBZfAN3P+HuHlvtBTgiItKpRAn8QcDrcev1sW0fY2bTzexF4HHgW4kOZGbzY1M+VWdOHmtPvSIi0k5RAt8SbPvEGby7P+buw4B/A+5OdCB3X+Xuo9x9VHbPvm0qVEREOiZK4NcDg+PW84CDLTV29+3ApWbWv4O1iYhIEkUJ/EpgqJnlm1l3YCZQEd/AzP6XmVlseSTQHTic7GJFRKT9Wr1Kx90bzWwh8ASQDax2971mtiC2fyXwVeB6M/sAaABmxH2JKyIinUCrgQ/g7luALc22rYxbXgosTW5pIiKSTLrTVkQkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEJEC38xKzazWzOrMrDzB/jlmVhP7edbMvpD8UkVEpCNaDXwzywYeAK4BCoBZZlbQrNkrwAR3LwTuBlYlu1AREemYbhHaFAN17v4ygJmtA8qAfR81cPdn49o/D+Qls0gRSY5Vq1axdu3adJchaRIl8AcBr8et1wNXfkr7/w38PtEOM5sPzAfI7vPZiCWKSLKsXbuW7c9X0X1AfrpLkTSIEviWYJsnbGj2JZoC/+pE+919FbHpnh4DhyY8hoikVvcB+Vw0e0m6y5B2OrB0WrtfGyXw64HBcet5wMHmjcysEPhP4Bp3P9zuikREJCWiXKVTCQw1s3wz6w7MBCriG5jZPwEbgX9395eSX6aIiHRUq2f47t5oZguBJ4BsYLW77zWzBbH9K4FFQD/gQTMDaHT3UakrW0RE2irKlA7uvgXY0mzbyrjlG4Ebk1uaiIgkk+60FREJhAJfRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAqHAFxEJhAJfRCQQCnwRkUAo8EVEAhEp8M2s1MxqzazOzMoT7B9mZs+Z2ftm9v3klykiIh3VrbUGZpYNPABMBuqBSjOrcPd9cc2OAN8F/i0VRYqISMe1GvhAMVDn7i8DmNk6oAw4G/jufgg4ZGb/mpIqpdOprq6mpKQk3WVIG1VXV0OfwekuQ9IkSuAPAl6PW68HrmzPm5nZfGA+QHafz7bnENIJ9CooYcSp3ekuQ9qhqKiI3bkj0l2GpEmUwLcE27w9b+buq4BVAD0GDm3XMST9eheV8iql6S5D2ql3uguQtInypW09EP87YB5wMDXliIhIqkQJ/EpgqJnlm1l3YCZQkdqyREQk2Vqd0nH3RjNbCDwBZAOr3X2vmS2I7V9pZhcBVUAf4EMzuxUocPd3U1e6iIi0RZQ5fNx9C7Cl2baVcctv0TTVIyIinZTutBURCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFAKPBFRAKhwBcRCYQCX0QkEAp8EZFAKPBFRAKhwBcRCUSkwDezUjOrNbM6MytPsN/MbEVsf42ZjUx+qSIi0hGtBr6ZZQMPANcABcAsMyto1uwaYGjsZz7wf5Jcp4iIdFCUM/xioM7dX3b308A6oKxZmzLg/3qT54ELzGxgkmsVEZEO6BahzSDg9bj1euDKCG0GAW/GNzKz+TT9BgDw/oGl0/a0qdrM0h94J91FpJD6l7m6ct+g6/fv8+19YZTAtwTbvB1tcPdVwCoAM6ty91ER3j8jqX+ZrSv3ryv3DcLoX3tfG2VKpx4YHLeeBxxsRxsREUmjKIFfCQw1s3wz6w7MBCqatakAro9drTMGOObubzY/kIiIpE+rUzru3mhmC4EngGxgtbvvNbMFsf0rgS3AVKAOOAnMjfDeq9pddWZQ/zJbV+5fV+4bqH8tMvdPTLWLiEgXpDttRUQCocAXEQlEygM/wmMZSszsmJlVx34WpbqmZDGz1WZ2yMwS3k+Q6Y+ciNC/TB67wWb2tJntN7O9ZnZLgjYZO34R+5fJ45drZjvN7C+x/v0kQZtMHr8o/Wv7+Ll7yn5o+pL3b8AlQHfgL0BBszYlwOZU1pHC/o0HRgJ7Wtg/Ffg9TfcpjAH+O901J7l/mTx2A4GRseXewEsJ/m1m7PhF7F8mj58B58eWc4D/BsZ0ofGL0r82j1+qz/CjPJYhY7n7duDIpzTJ6EdOROhfxnL3N939hdjycWA/TXeHx8vY8YvYv4wVG5MTsdWc2E/zK1Ayefyi9K/NUh34LT1yobmxsV9dfm9mV6S4pnMpav8zWcaPnZkNAf6ZprOoeF1i/D6lf5DB42dm2WZWDRwC/uDuXWr8IvQP2jh+qQ78KI9ceAG42N2/ANwHbEpxTedSpEdOZLCMHzszOx94FLjV3d9tvjvBSzJq/FrpX0aPn7ufcfcimu7sLzaz4c2aZPT4Rehfm8cv1YHf6iMX3P3dj351cfctQI6Z9U9xXedKl37kRKaPnZnl0BSGj7j7xgRNMnr8Wutfpo/fR9z9KLANKG22K6PH7yMt9a8945fqwG/1sQxmdpGZWWy5OFbT4RTXda506UdOZPLYxer+FbDf3X/WQrOMHb8o/cvw8fusmV0QWz4PmAS82KxZJo9fq/1rz/hFeVpmu3m0xzJ8DbjJzBqBBmCmx76C7uzM7Lc0fVPe38zqgR/R9OXKR31rzyMnOo0I/cvYsQOuAv4d2B2bJwW4A/gn6BLjF6V/mTx+A4GHrekPNGUB/+Xum63jj3zpLKL0r83jp0criIgEQnfaiogEQoEvIhIIBb6ISCAU+CIigVDgi4gEQoEvIhIIBb6ISCD+B3rpvNSBZOBSAAAAAElFTkSuQmCC\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/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAc8UlEQVR4nO3df3BU9f3v8ecbBFMpqEVpO4ASWi2NJN3aEAFb2RawqVW5TuuAsa21VyNt1TpTvVK/lrlfy1dxdPrDH/fSndbqbd3Kd/xVxlL1azVQrRSC3UaIpaWAJVoav1gxIDEi7/vHJsvJsiGbZDdnc/J6zDjm7DnZfXsM73x4n/fn8zF3R0REhr4RYQcgIiKFoYQuIhIRSugiIhGhhC4iEhFK6CIiEXFUWB98wgkn+JQpU8L6eBGRIWnjxo3/7e4n5joXWkKfMmUKjY2NYX28iMiQZGYv93ROJRcRkYhQQhcRiQgldBGRiAithi4i0fDOO+/Q0tJCe3t72KFESllZGZMmTWLUqFF5f48SuogMSEtLC2PHjmXKlCmYWdjhRIK7s3v3blpaWigvL8/7+1RyEZEBaW9vZ/z48UrmBWRmjB8/vs9/61FCF5EBUzIvvP7cUyV0EZGIUEIXEemHVCrF6tWr+/x98Xi8aJMq9VBUBJiy5Nfdjncs/3xIkchQkUqlaGxs5Jxzzgk7lAyN0EUkEn7xi19QU1NDLBbjiiuu4A9/+ANVVVW0t7ezb98+TjvtNDZt2kRDQwNnnXUWF1xwARUVFSxevJiDBw8C8OSTTzJr1ixOP/10LrzwQvbu3QvAhg0bmD17Nh/72Meoqalhz549LF26lJUrVxKLxVi5ciX79u3ja1/7GjNmzODjH/84v/rVrwDYv38/ixYtoqqqioULF7J///6i3QON0EWkYK655hpSqVRB3zMWi/HDH/7wiNe89NJLrFy5kueee45Ro0bxjW98gy1btnD++edz4403sn//fr70pS8xffp0GhoaWL9+Pc3NzZx88snU1tby8MMPE4/HWbZsGU899RRjxozh1ltv5fvf/z5Llixh4cKFrFy5khkzZvDmm29yzDHHcNNNN9HY2Mhdd90FwA033MBnPvMZ7rnnHt544w1qamqYN28eP/7xjznmmGNoamqiqamJ008/vaD3J0gJXUSGvN/+9rds3LiRGTNmAOlR8YQJE1i6dCkzZsygrKyMO+64I3N9TU0NU6dOBeCiiy7i2WefpaysjObmZs4880wAOjo6mDVrFlu2bOGDH/xg5r3HjRuXM4Ynn3ySVatWcfvttwPpds6///3vrF27lquvvhqAqqoqqqqqinMTUEIXkQLqbSRdLO7OJZdcwi233NLt9V27drF3717eeecd2tvbGTNmDHB4S6CZ4e7Mnz+fX/7yl93ONTU15dVC6O489NBDfOQjHzns3GC1deZVQzezWjPbYmZbzWxJjvNxM9tjZqnOf5YWPlQRkdzmzp3Lgw8+SGtrKwCvv/46L7/8MvX19Xzve9/j4osv5vrrr89cv379erZv387BgwdZuXIln/zkJ5k5cybPPfccW7duBeCtt97iL3/5C9OmTePVV19lw4YNALS1tXHgwAHGjh1LW1tb5j0/+9nPcuedd+LuAPzxj38E4KyzzuL+++8HYNOmTTQ1NRXtPvQ6QjezkcDdwHygBdhgZqvcvTnr0t+5+7lFiFFE5IgqKipYtmwZZ599NgcPHmTUqFEsWLCAo446irq6Ot59911mz57N008/zYgRI5g1axZLlizhxRdfzDwgHTFiBPfeey8XXXQRb7/9NgDLli3j1FNPZeXKlVx11VXs37+f97znPTz11FN8+tOfZvny5cRiMb7zne/w3e9+l2uuuYaqqircnSlTpvDYY4/x9a9/nUsvvZSqqipisRg1NTVFuw/W9dukxwvMZgH/290/23n8HQB3vyVwTRy4ti8Jvbq62rXBhZQKtS3230svvcRHP/rRsMPIW0NDA7fffjuPPfZY2KH0Kte9NbON7l6d6/p8Si4TgZ2B45bO17LNMrM/mdlvzOy0XG9kZvVm1mhmja+99loeHy0iIvnKJ6HnquZnD+tfAE52948BdwKP5nojd0+4e7W7V594Ys4t8UREiioejw+J0Xl/5JPQW4DJgeNJwKvBC9z9TXff2/n1amCUmZ1QsChFRKRX+ST0DcApZlZuZqOBRcCq4AVm9gHr7Msxs5rO991d6GBFRKRnvXa5uPsBM7sSeAIYCdzj7pvNbHHn+RXAF4Gvm9kBYD+wyHt72ioiIgWV18SizjLK6qzXVgS+vgu4q7ChiYhIX2imqIgUVHYL6EDl00K6Y8cOzj33XDZt2tTt9aVLl3LWWWcxb968nN/36KOPcuqpp1JRUVGQWMOm1RZFJLJuuummHpM5pBN6c3P2HMn+OXDgQEHeZyCU0EUkEt59910uv/xyTjvtNM4++2z279/PV7/6VR588EEAlixZQkVFBVVVVVx77bX8/ve/Z9WqVVx33XXEYjH+9re/kUqlmDlzJlVVVVxwwQX861//AtLL51ZVVTFr1iyuu+46pk+fDsC9997LhRdeyHnnncfZZ5/N3r17mTt3LqeffjqVlZWZJXR37NjBtGnTuOyyy5g+fToXX3wxTz31FGeeeSannHIK69evL8g9UEIXkUj461//yje/+U02b97Mcccdx0MPPZQ59/rrr/PII4+wefNmmpqauPHGG5k9ezbnn38+t912G6lUig996EN85Stf4dZbb6WpqYnKykr+/d//HYBLL72UFStW8PzzzzNy5Mhun/v8889z33338fTTT1NWVsYjjzzCCy+8wDPPPMO3v/3tzNouW7du5Vvf+hZNTU38+c9/JplM8uyzz3L77bdz8803F+QeKKGLSCSUl5cTi8UA+MQnPsGOHTsy58aNG0dZWRmXXXYZDz/8MMccc8xh379nzx7eeOMN5syZA8All1zC2rVreeONN2hra2P27NkA1NXVdfu++fPn8773vQ9Ir7h4ww03UFVVxbx583jllVf45z//mYmvsrKSESNGcNpppzF37lzMjMrKym6xDoQSuohEwtFHH535euTIkd1q2kcddRTr16/nC1/4Ao8++ii1tbV5v29vHdhdS/IC3H///bz22mts3LiRVCrF+9//ftrb2w+Lb8SIEZnjESNGFKz+ri4XiZxgl8WROiQK3Y0hpWvv3r289dZbnHPOOcycOZMPf/jDAN2WwD322GM5/vjj+d3vfsenPvUpfv7znzNnzhyOP/54xo4dy7p165g5cyYPPPBAj5+zZ88eJkyYwKhRo3jmmWd4+eWXB+W/r4sSuogUVCmuVNnW1saCBQtob2/H3fnBD34AwKJFi7j88su54447ePDBB7nvvvtYvHgxb731FlOnTuVnP/sZAD/96U+5/PLLGTNmDPF4nGOPPTbn51x88cWcd955VFdXE4vFmDZt2qD9N0Iey+cWi5bPlWIpxAi9FJNSqRpqy+f2x969e3nve98LwPLly/nHP/7Bj370o6J/bl+Xz9UIXUSkF7/+9a+55ZZbOHDgACeffDL33ntv2CHlpIQuItKLhQsXsnDhwrDD6JW6XERkwLQWX+H1554qoYvIgJSVlbF7924l9QJyd3bv3k1ZWVmfvk8lFxEZkEmTJtHS0oK2lSyssrIyJk2a1KfvUUIXkQEZNWoU5eXlYYchqOQiIhIZSugiIhGhhC4iEhGqoUskFHNdlnxnnoqETSN0EZGIUEIXEYkIJXQRkYhQQhcRiQgldBGRiFBCFxGJCCV0EZGIUB+6RFp2f7r6yCXKNEIXEYkIJXQRkYhQQhcRiYi8ErqZ1ZrZFjPbamZLjnDdDDN718y+WLgQRUQkH70mdDMbCdwNfA6oAC4ys4oerrsVeKLQQYqISO/yGaHXAFvdfZu7dwAPAAtyXHcV8BDQWsD4REQkT/m0LU4EdgaOW4AzgheY2UTgAuAzwIyCRScyCNpSj7OvuaHba/F1tx06X1bJ2FjtIEcl0nf5jNAtx2vZ23v/ELje3d894huZ1ZtZo5k1akNZKRX7mhvoaN2e81wqlTos2YuUqnxG6C3A5MDxJODVrGuqgQfMDOAE4BwzO+DujwYvcvcEkACorq7O/qUgMmgSiQS7kncC0NG6ndETyvlA3fLM+YbOCUjxeJy16xrZlUz3AiSmvkJ9ff3gByySh3xG6BuAU8ys3MxGA4uAVcEL3L3c3ae4+xTgQeAb2clcpJQkk8nMqHz0hHLGVMRzXldXV8foCekd7Ttat5NMJgcrRJE+63WE7u4HzOxK0t0rI4F73H2zmS3uPL+iyDGKFEX2qDyX+vp6bt42ESAzShcpVXmt5eLuq4HVWa/lTOTu/tWBhyVSeG2px4nH0w87U6kUjJt85G/IIZVKEY/HM8d1dXUqwUjJ0ExRGTb2NTekEzkQi8V6LLP0ZExFnFgsljlOpVIqwUhJ0WqLMqzEYjEaGhqAw1di7M3YWC0Ny+/MHAdH6iKlQCN0EZGI0AhdIit7wlBH63aYOr6gnxGsqaueLmFTQpfI6pow1NV2OHpCOXV1dQV7/+B7ddXmldAlTEroEmnZrYn19YXbsai+vj6TwFVPl1KghC7SB9rSTkqZHoqKiESEErqISESo5CKREuxsCT4QHQzqeJGwKaFLpAQ7W3ItutXXyUT5UseLlAIldImcfBbdKjR1vEgpUA1dRCQilNBFRCJCJReRHIpVaxcpJo3QRUQiQiN0GfJy7Q8aNm2EIWFQQpchr2t/0J5aFQdb9gJgamOUwaKELpEQRqtitkN194ns6NxEA9TGKINHNXQRkYhQQhcRiQgldBGRiFBCFxkE67btZsqSX6u/XYpKCV1EJCLU5SJDTiKRIJlMZo5TqRSMmxxeQHnoaN3OruQSABJTX1ELoxSFRugy5CSTyUxvN0AsFgu99/xI6urqMpOdOlq3d/tlJFJIGqHLkBSLxWgI9HqXcm26vr6em7dNBMiM0kWKQSN0EZGIUEIXEYkIJXQRkYhQDV1kAEq5di/DT14jdDOrNbMtZrbVzA57qmNmC8ysycxSZtZoZp8sfKgiInIkvY7QzWwkcDcwH2gBNpjZKndvDlz2W2CVu7uZVQH/CUwrRsAiIpJbPiP0GmCru29z9w7gAWBB8AJ33+vu3nk4BnBECiiRSBCPx4nH49160Ieirs0v4vE4iUQi7HAkQvKpoU8EdgaOW4Azsi8yswuAW4AJwOdzvZGZ1QP1ACeddFJfY5VhLJlMsnZdY3qCzrjJh20iMVSMqYhT2f4ioI0vpPDySeiW47XDRuDu/gjwiJmdBXwPmJfjmgSQAKiurtYoXvokuIlFfX3OMUPJGxurpWF5ers8bXwhhZZPyaUFCC6UMQl4taeL3X0t8CEzO2GAsYmISB/kk9A3AKeYWbmZjQYWAauCF5jZh83MOr8+HRgN7C50sCIi0rNeSy7ufsDMrgSeAEYC97j7ZjNb3Hl+BfAF4Ctm9g6wH1gYeEgqUnDq/xY5XF4Ti9x9NbA667UVga9vBW4tbGgiItIXmvovIhIRSugiIhGhtVykpHTVxttSj2f6tWFo7EoUlG+Nf+26RspOqswc33HjVepLl37TCF1K0r7mhiG1K1F/BHcyAu1mJAOnEbqUrKG0K1F/BHcyAu1mJAOnEbqISEQooYuIRIQSuohIRCihi4hEhBK6iEhEKKGLiESEErqISESoD11KRiKRYFcyvflDR+t2mDo+5IgGX9f2dJCeeKRZo9IXGqFLyUgmk+lETnp3oqG6zVx/jamIE4vFgHRi16xR6SuN0KWkBLeZu3kb3Byx2aFHou3pZKA0QhcRiQgldBGRiFBCFxGJCNXQZdBlr5q4Y/nnQ4pEJFo0QhcRiQgldBGRiFDJRULVlnqcePw2YOhtM1dsmmQkfaWELqHa19zA2tbt6a3Yxk2O3DZzueSz81JwUlXXVnxK6NIbJXQJXXAykaTV19dnErgmGUm+VEMXEYkIJXQRkYhQQhcRiQgldBGRiFBCFxGJCCV0EZGIUNuiyBAU7GXXWjjSJa+Ebma1wI+AkcBP3H151vmLges7D/cCX3f3PxUyUImG4DZzkN5qbvSE8hAjGhqCs0YB2soqGRurDS8gKUm9JnQzGwncDcwHWoANZrbK3ZsDl20H5rj7v8zsc0ACOKMYAcvQ1rXNXFcSHz2hfFjMDh2I7K34UqkU7eN2K6HLYfIZodcAW919G4CZPQAsADIJ3d1/H7h+HTCpkEFKtGhmaN8EZ41Ceuboum27Q4xISlU+D0UnAjsDxy2dr/XkfwK/yXXCzOrNrNHMGl977bX8oxQRkV7lk9Atx2ue80KzT5NO6NfnOu/uCXevdvfqE088Mf8oRUSkV/mUXFqA4Jqmk4BXsy8ysyrgJ8Dn3F1/HxQRGWT5jNA3AKeYWbmZjQYWAauCF5jZScDDwJfd/S+FD1NERHrT6wjd3Q+Y2ZXAE6TbFu9x981mtrjz/ApgKTAe+D9mBnDA3auLF7YMNV1907v0MK/g1JMuXfLqQ3f31cDqrNdWBL6+DLissKGJiEhfaKaoFF1wMpEmEhVGR+t2diWXADCmIq6edAG0losMgq7JRKCJRIVQV1eX+aXY0bqdfc0N4QYkJUMjdCmY7L0yg/VcTSbKT0/7jQbvZX19PTdvS08F6Rqli4BG6CIikaGELiISEUroIiIRoRq6FI16z0UGl0boIiIRoYQuIhIRKrlIUbSlHs/0R2syUXEFJxklpr7Sbe10GV6U0KUo9jU3ZBK5JhMVT/C+drRu5+pld2Z61EFruww3SuhSNJpMVHxjY7WZaf+aZCSqoYuIRIQSuohIRCihi4hEhBK6iEhEKKGLiESEulykIIKbWIB6z8MS7EkH9aUPNxqhS0EEN7EAbWQRhjEV8W6/RDtat5NMJkOMSAabRuhSMOo7D1ewJx3Ulz4caYQuIhIRSugiIhGhkovIENDTXqMiQRqhi4hEhBK6iEhEKKFLvyUSCeLxOPF4nFQqFXY4kkMqlcr8P0okEmGHI0WmGrr0WzKZZO26xnTv87jJ6jsvMWMq4uxrbmDdtt2ZOQLBSUbBurzWTY8GJXQZEPWely6tlT78qOQiIhIRSugiIhGhhC4iEhF5JXQzqzWzLWa21cwOK8aZ2TQze97M3jazawsfpoiI9KbXh6JmNhK4G5gPtAAbzGyVuzcHLnsduBr4H8UIUkpHIpHIrOCXSqVg3ORwA5K8dbUwdmkrq+y2mJcMffmM0GuAre6+zd07gAeABcEL3L3V3TcA7xQhRikhyWQy03Mei8XUqjhEjKmIE4vFMsepVIp9zQ2hxSPFkU/b4kRgZ+C4BTijPx9mZvVAPcBJJ53Un7eQEtA+bjI7Zl4HwNiQY5H8jI3V0rD80AYk8Xicddt2hxiRFEM+I3TL8Zr358PcPeHu1e5efeKJJ/bnLUREpAf5JPQWIFgonQS8WpxwRESkv/JJ6BuAU8ys3MxGA4uAVcUNS0RE+qrXGrq7HzCzK4EngJHAPe6+2cwWd55fYWYfABqBccBBM7sGqHD3N4sXugyGYFcLqLMlSoIbSmsz6WjIqw/d3Ve7+6nu/iF3/4/O11a4+4rOr3e5+yR3H+fux3V+rWQeAcGuFlBnS1TU1dVlNpTWZtLRocW5pFexWIyGhobMsXbPGfrq6+u5edtEQAt3RYmm/ouIRIRG6JKh9bGjTX+zij6N0EVEIkIJXQ7TlnpcW8sNM9qqLhqU0OUw+5obuq3XUldXF25AUlTBdV5SqZQ6XoYw1dCHsSPVVLM7WyS6guu8BFdjlKFHI3QRkYhQQhcRiQiVXARIPwjtWh+7o3U7TB0fbkASmuyNMOrq6rQswBChhC5A+kFoR+t2Rk8oZ/SE8m4PQtW/PHxkPwDvejieK6Fn/1xo7kL4lNAlY/SEcj5QtxyA+nr94RyO6uvruyVvPSQdWlRDFxGJCI3Qh6lEIsGu5KEtybrKLSLZgjV11dNLmxL6MNNV99yVvLNbEh89oVzL4sph6urqWLdtN+u27U4/LCd3PV1KgxL6MBasmYvkomV2hxbV0EVEIkIj9GEkWDdXzVz6I1hPbyurZGysNtyApBsl9GEkmUx26zVXzVy65DPXYExFnMr2F4F0Ym8ft1sJvcQooQ8zqptLf2Uv4rVu2+6QI5JsSugRl0gkMsuhplIpGDc53IAkMjpat3d7UJqY+oo6YEKmh6IRl0wmu61trjKLFEJdXV23ZzAdrdu1jnoJ0Ah9GAiuba51WaQQgu2MUJiWRq0NM3BK6BETLLFAuszStRuNSDFpRmn4VHKJmGCJBbSFnAwObWNXGjRCj4DsB5/aPk4GW3YHjEbr4VBCH0J6qjF2jcpjsRixWIwXyyq7XatapBTKkZ7BdJ1rK6ukq8qXaz11/WwWjxL6ENWWepx4/Dbg8FG5HnxKmI40WgfNMC0mJfQhInu527d3bmINMGfOHNXJpWRl/1yuWbMGWJPZ7lC964WlhF7CgrXx9B8EOHry9My/77jxKv1hkJKWvQNSIpHg6mXpgcnbOzdxxRVXZH7GNXIfOCX0QZJvGSS4WfPbOzcB6VH4nDlzeDHrBz64TVyhyywq20gxBPvXu37W123b3fmzfmjkDjA+9XjOBK+6e8/ySuhmVgv8CBgJ/MTdl2edt87z5wBvAV919xcKHGtkBJN2tq4kfvTk6YeNwpVkJUrGxmozCTv7z8TbOzfx9s5NOf+cqEzTs14TupmNBO4G5gMtwAYzW+XuzYHLPgec0vnPGcD/7fx3pGVP4jmSXYGFjIJJO9vRk6czpiKe+UHXZs0yHASTO/Q86Mku0/TFcGifzGeEXgNsdfdtAGb2ALAACCb0BcD/c3cH1pnZcWb2QXf/R09vumXLliG/o3hXXXvOnDl9+r7spC0i3WUn+C5tqcczS/j2xZo1a1izZk3kJzxZOgcf4QKzLwK17n5Z5/GXgTPc/crANY8By9392c7j3wLXu3tj1nvVA12/Ij8CbCnUf8gAnAD8d9hBlAjdi0N0Lw7RvTikFO7Fye5+Yq4T+YzQLcdr2b8F8rkGd08AiTw+c9CYWaO7V4cdRynQvThE9+IQ3YtDSv1e5LOWSwsQXER7EvBqP64REZEiyiehbwBOMbNyMxsNLAJWZV2zCviKpc0E9hypfi4iIoXXa8nF3Q+Y2ZXAE6TbFu9x981mtrjz/ApgNemWxa2k2xYvLV7IBVdSJaCQ6V4contxiO7FISV9L3p9KCoiIkOD1kMXEYkIJXQRkYhQQu9kZteamZvZCWHHEhYzu83M/mxmTWb2iJkdF3ZMg83Mas1si5ltNbOBb5Q5RJnZZDN7xsxeMrPNZvatsGMKm5mNNLM/ds67KUlK6KR/eEkvbfD3sGMJ2X8B0929CvgL8J2Q4xlUgWUuPgdUABeZWUW4UYXmAPBtd/8oMBP45jC+F12+BbwUdhBHooSe9gPgf5FjMtRw4u5PuvuBzsN1pOcTDCeZZS7cvQPoWuZi2HH3f3QtsOfubaQT2cRwowqPmU0CPg/8JOxYjmTYJ3QzOx94xd3/FHYsJeZrwG/CDmKQTQR2Bo5bGMZJrIuZTQE+Dvwh5FDC9EPSg76DIcdxRMNiPXQzewr4QI5T/wbcAJw9uBGF50j3wt1/1XnNv5H+K/f9gxlbCchrCYvhxMzeCzwEXOPub4YdTxjM7Fyg1d03mlk85HCOaFgkdHefl+t1M6sEyoE/pZd0ZxLwgpnVuPuuQQxx0PR0L7qY2SXAucBcH36TFLSERYCZjSKdzO9394fDjidEZwLnm9k5QBkwzsx+4e5fCjmuw2hiUYCZ7QCq3T3s1dRC0bmRyfeBOe7+WtjxDDYzO4r0w+C5wCukl72oc/fNoQYWgs5Na+4DXnf3a0IOp2R0jtCvdfdzQw4lp2FfQ5du7gLGAv9lZikzWxF2QIOp84Fw1zIXLwH/ORyTeaczgS8Dn+n8WUh1jlClhGmELiISERqhi4hEhBK6iEhEKKGLiESEErqISEQooYuIRIQSuohIRCihi4hExP8H0iUfIFEdSf8AAAAASUVORK5CYII=\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"def sample_Zn(alpha,b,n):\n",
|
|
" assert alpha > 2\n",
|
|
" assert b > 0\n",
|
|
" assert n >= 1 and type(n) == int\n",
|
|
" \n",
|
|
" E_X = alpha*b/(alpha - 1)\n",
|
|
" Var_X = alpha*b**2/( (alpha - 1)**2*(alpha - 2) )\n",
|
|
" \n",
|
|
" inv_pareto_samples = [inversion_sample(lambda p: f_inv_pareto(alpha, b, p)) for _ in range(n)]\n",
|
|
" return np.sqrt(n/Var_X)*(np.mean(inv_pareto_samples) - E_X)\n",
|
|
"\n",
|
|
"# Plotting\n",
|
|
"alpha = 4\n",
|
|
"b = 1\n",
|
|
"n = 1000\n",
|
|
"pdf = gaussian\n",
|
|
"samples = [inversion_sample(lambda p: sample_Zn(alpha, b, n)) for _ in range(1000)]\n",
|
|
"compare_plot(samples, pdf, -5, 5, 100)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 10,
|
|
"id": "b5360d77",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "e50b33644ddd6bce391b36cefcc2e308",
|
|
"grade": true,
|
|
"grade_id": "cell-5d16b014bef9d86f",
|
|
"locked": true,
|
|
"points": 10,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"assert_almost_equal(np.mean([sample_Zn(3.5,2.1,100) for _ in range(100)]),0,delta=0.3)\n",
|
|
"assert_almost_equal(np.std([sample_Zn(3.5,2.1,100) for _ in range(100)]),1,delta=0.3)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "6192f05d",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "08ece68d59de21d798d9a955f59be690",
|
|
"grade": false,
|
|
"grade_id": "cell-3e7a23657e9b8374",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"__(b)__ Now take $\\alpha = 3/2$ and $b=1$. \n",
|
|
"With some work (which you do not have to do) one can show that the characteristic function of $X$ admits the following expansion around $t=0$,\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\varphi_X(t) = 1 + 3 i t - (|t|+i t)\\,\\sqrt{2\\pi|t|} + O(t^{2}).\n",
|
|
"$$\n",
|
|
"\n",
|
|
"Based on this, prove the **generalized CLT** for this particular distribution $X$ which states that $Z_n = c\\, n^{1/3} (\\bar{X}_n - \\mathbb{E}[X])$ in the limit $n\\rightarrow\\infty$ converges in distribution, with a to-be-determined choice of overall constant $c$, to a limiting random variable $\\mathcal{S}$ with characteristic function \n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\varphi_{\\mathcal{S}}(t) = \\exp\\big(-(|t|+it)\\sqrt{|t|}\\big).\n",
|
|
"$$\n",
|
|
"\n",
|
|
"__(15 pts)__"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "9735cd88",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "dfd8683eea5663baa81f138a2809722b",
|
|
"grade": true,
|
|
"grade_id": "cell-b25551eca32c4807",
|
|
"locked": false,
|
|
"points": 15,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"\\begin{align}\n",
|
|
"\\phi_{Z_n}(t) &= \\mathbb{E}\\left[ e^{itZ_n} \\right]\n",
|
|
" \\\\ &= \\mathbb{E}\\left[ e^{itcn^{1/3}(\\bar{X_n} - \\mathbb{E}[X])} \\right]\n",
|
|
" \\\\ &= \\mathbb{E}\\left[ e^{itcn^{1/3}(\\frac{1}{n}\\sum_{i=1}^n X_i - \\mathbb{E}[X])} \\right]\n",
|
|
" \\\\ &= \\mathbb{E}\\left[ \\left( \\prod_{i=1}^n e^{itcn^{-2/3}X_i} \\right) e^{itcn^{1/3}\\mathbb{E}[X])} \\right]\n",
|
|
" \\\\ &= \\left( \\prod_{i=1}^n \\mathbb{E}\\left[ e^{itcn^{-2/3}X_i} \\right] \\right)\\mathbb{E}\\left[ e^{itcn^{1/3}\\mathbb{E}[X])} \\right]\n",
|
|
" \\\\ &= \\left( \\prod_{i=1}^n \\phi_X(cn^{-2/3}t) \\right)\\mathbb{E}\\left[ e^{itcn^{1/3}\\mathbb{E}[X])} \\right]\n",
|
|
" \\\\ &= \\left( \\phi_X(cn^{-2/3}t) \\right)^n \\mathbb{E}\\left[ e^{itcn^{1/3}\\mathbb{E}[X])} \\right]\n",
|
|
"\\end{align}\n",
|
|
"where we used the identity for products of indepedent expectation values <https://hef.ru.nl/~tbudd/mct/lectures/probability_random_variables.html#equation-product-expectation>, and the definition of $\\phi_X(t) := \\mathbb{E}\\left[ e^{itX} \\right]$.\n",
|
|
"\n",
|
|
"Next, we will use the Taylor expansion around $t = 0$ as is given above, and, for the latter exponential, $\\mathbb{E}(X) = 3$ for $\\alpha = 3/2, b = 1$ as given.\n",
|
|
"\n",
|
|
"\\begin{align}\n",
|
|
"\\phi_{Z_n}(t) &= \\left( \\phi_X(cn^{-2/3}t) \\right)^n \\mathbb{E}\\left[ e^{itcn^{1/3}\\mathbb{E}[X])} \\right]\n",
|
|
" \\\\ &= \\left( 1 + 3 i cn^{-2/3}t - (|cn^{-2/3}t|+i cn^{-2/3}t)\\,\\sqrt{2\\pi|cn^{-2/3}t|} + \\mathcal{O}(t^2) \\right)^n e^{3itcn^{1/3}}\n",
|
|
" \\\\ &= \\left( 1 + \\frac{1}{n} \\left[ 3 i cn^{1/3}t - (|ct|+i ct)\\,\\sqrt{2\\pi|ct|} \\right] + \\mathcal{O}(t^2) \\right)^n e^{3itcn^{1/3}}\n",
|
|
"\\end{align}\n",
|
|
"\n",
|
|
"Taking the limit $n \\to \\infty$, the first set of parentheses can be rewritten in terms of an exponential using the identity $\\lim_{n\\to\\infty} (1 + \\frac{a}{n})^n = e^{a}$, we find a way to our desired expression.\n",
|
|
"\n",
|
|
"\\begin{align}\n",
|
|
"\\lim_{n\\to\\infty} \\phi_{Z_n}(t) &= \\lim_{n\\to\\infty} \\left( 1 + \\frac{1}{n} \\left[ 3 i cn^{1/3}t - (|ct|+i ct)\\,\\sqrt{2\\pi|ct|} \\right] + \\mathcal{O}(t^2) \\right)^n e^{-3itcn^{1/3}}\n",
|
|
" \\\\ &= \\lim_{n\\to\\infty} \\exp{({3 i cn^{1/3}t - (|ct|+i ct)\\,\\sqrt{2\\pi|ct|}})} \\exp{(e^{-3itcn^{1/3}})}\n",
|
|
" \\\\ &= \\exp{({-(|ct|+i ct)\\,\\sqrt{2\\pi|ct|}})}\n",
|
|
"\\end{align}\n",
|
|
"\n",
|
|
"This matches $\\phi_S(t) = \\exp\\big(-(|t|+it)\\sqrt{|t|}\\big)$ for $\\sqrt{2\\pi} c^{3/2} = 1 \\implies c = (2\\pi)^{\\frac{-1}{3}}$."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "5b1d9f54",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"editable": false,
|
|
"nbgrader": {
|
|
"cell_type": "markdown",
|
|
"checksum": "342020128f929d47eabfdf9c075ff20c",
|
|
"grade": false,
|
|
"grade_id": "cell-d1701433c3c77172",
|
|
"locked": true,
|
|
"schema_version": 3,
|
|
"solution": false,
|
|
"task": false
|
|
}
|
|
},
|
|
"source": [
|
|
"__(c)__ The random variable $\\mathcal{S}$ has a [stable Lévy distribution](https://en.wikipedia.org/wiki/Stable_distribution) with index $\\alpha = 3/2$ and skewness $\\beta = 1$. Its probability density function $f_{\\mathcal{S}}(x)$ does not admit a simple expression, but can be accessed numerically using SciPy's `scipy.stats.levy_stable.pdf(x,1.5,1.0)`. Verify numerically that the generalized CLT of part (b) holds by comparing an appropriate histogram to this PDF. __(10 pts)__"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 11,
|
|
"id": "b06896e5",
|
|
"metadata": {
|
|
"deletable": false,
|
|
"nbgrader": {
|
|
"cell_type": "code",
|
|
"checksum": "c6fe081427f342c354ee8a9b3b3331e7",
|
|
"grade": true,
|
|
"grade_id": "cell-e08d054985cfa762",
|
|
"locked": false,
|
|
"points": 10,
|
|
"schema_version": 3,
|
|
"solution": true,
|
|
"task": false
|
|
}
|
|
},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD7CAYAAABkO19ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAgi0lEQVR4nO3de3BU9f3/8eebW6NpvBTBdkQEWyyNJK40IJeOhCo0XqnTrwPGVmu/Gmi91E71V+rXOlPLWJw6rUr9fjFjrX5bU/mOimWUqrUa0CrlYtMIKC0FWiNSKCoGJSL6/v2xm+Vk2SQnyW529+zrMcOwe875bD674nvfeZ/PxdwdERGJrgG57oCIiGSXAr2ISMQp0IuIRJwCvYhIxCnQi4hEnAK9iEjEhQr0ZlZjZpvMbLOZzU9zfpaZNZtZk5mtNbMvBM5tM7OX289lsvMiItI9624cvZkNBP4KzABagDXARe6+MXDNx4F33d3NrBL4P3cfmzi3Dahy939n5y2IiEhXBoW4ZiKw2d23AJjZg8AsIBno3X1v4PpSoE+zsI455hgfNWpUX15CRKSorFu37t/uPizduTCB/jjgtcDzFuC01IvM7ALgx8Bw4JzAKQeeMjMH7nb3+u5+4KhRo1i7VlUeEZGwzOwfnZ0LU6O3NMcOydjdfWmiXPNl4EeBU1PdfTxwFnClmZ3eSSfrEvX9tbt27QrRLRERCSNMoG8Bjg88HwFs7+xid18JfNrMjkk83574eyewlHgpKF27enevcveqYcPS/vYhIiK9ECbQrwHGmNloMxsCzAGWBS8ws8+YmSUejweGALvNrNTMyhLHS4GZwPpMvgEREelatzV6dz9gZlcBTwIDgXvdfYOZzUucXwx8BbjEzD4A9gGzEyNwjgWWJr4DBgEN7v5Elt6LiOSZDz74gJaWFtra2nLdlcgoKSlhxIgRDB48OHSbbodX5kJVVZXrZqxI4du6dStlZWUMHTqURMInfeDu7N69m9bWVkaPHt3hnJmtc/eqdO00M1ZEsqatrU1BPoPMjKFDh/b4NyQFehHJKgX5zOrN56lALyKSQU1NTSxfvrzH7aqrq7M2f0iBXkQkg3ob6LNJgV7yyqj5jyf/iGTCr3/9ayZOnEgsFmPu3Ln86U9/orKykra2Nt59911OPvlk1q9fT2NjI6effjoXXHAB5eXlzJs3j48++giAp556ismTJzN+/HguvPBC9u6Nr/qyZs0apkyZwimnnMLEiRPZs2cPN910E0uWLCEWi7FkyRLeffddvvGNbzBhwgROPfVUfvvb3wKwb98+5syZQ2VlJbNnz2bfvn1Z+wzCLIEgIlKQXnnlFZYsWcIf//hHBg8ezLe+9S02bdrE+eefz4033si+ffv46le/yrhx42hsbGT16tVs3LiRE044gZqaGh555BGqq6tZsGABTz/9NKWlpdx666389Kc/Zf78+cyePZslS5YwYcIE3nnnHQ4//HBuvvlm1q5dy89//nMAbrjhBr74xS9y77338vbbbzNx4kTOPPNM7r77bg4//HCam5tpbm5m/PjxWfscFOhFpF9ce+21NDU1ZfQ1Y7EYt99+e6fn//CHP7Bu3TomTJgAxLPo4cOHc9NNNzFhwgRKSkq48847k9dPnDiRE088EYCLLrqI559/npKSEjZu3MjUqVMB2L9/P5MnT2bTpk186lOfSr72EUcckbYPTz31FMuWLeO2224D4iOR/vnPf7Jy5UquueYaACorK6msrOzbh9EFBXoRiSx359JLL+XHP/5xh+M7duxg7969fPDBB7S1tVFaWgocOqLFzHB3ZsyYwW9+85sO55qbm0ONgHF3Hn74YT772c8ecq6/RiQp0ItIv+gq886WM844g1mzZvGd73yH4cOH8+abb9La2srVV1/Nj370I7Zu3cr3vve9ZJll9erVbN26lRNOOIElS5ZQV1fHpEmTuPLKK9m8eTOf+cxneO+992hpaWHs2LFs376dNWvWMGHCBFpbWznssMMoKyujtbU12YcvfelLLFq0iEWLFmFm/PnPf+bUU0/l9NNP54EHHmD69OmsX7+e5ubmrH0OCvSSU7rpKtlUXl7OggULmDlzJh999BGDBw9m1qxZDBo0iNraWj788EOmTJnCM888w4ABA5g8eTLz58/n5ZdfTt6YHTBgAPfddx8XXXQR77//PgALFizgpJNOYsmSJVx99dXs27ePww47jKeffprp06ezcOFCYrEY3//+9/nBD37AtddeS2VlJe7OqFGjeOyxx/jmN7/JZZddRmVlJbFYjIkT0673mBFaAkFyqqtAv23hOZ2ek8Lwyiuv8LnPfS7X3QilsbGR2267jcceeyzXXelWus9VSyCIiBQxlW5ERIjPTK2urs51N7JCGb2ISMQpo5d+EazFq/Yu0r+U0YuIRJwCvYhIxKl0I/2uN2PnU9uo/CNhbdu2jXPPPZf16ztuV33TTTdx+umnc+aZZ6Zt9+ijj3LSSSdRXl7eH93MKgV6Eek3mZ4g15cv/JtvvrnL848++ijnnntuRgL9gQMHGDQod+FWpRvJW1qyWDLlww8/5IorruDkk09m5syZ7Nu3j69//es89NBDAMyfP5/y8nIqKyu57rrreOGFF1i2bBnXX389sViMv//97zQ1NTFp0iQqKyu54IILeOutt4D4UsWVlZVMnjyZ66+/nnHjxgFw3333ceGFF3Leeecxc+ZM9u7dyxlnnMH48eOpqKhILle8bds2xo4dy+WXX864ceO4+OKLefrpp5k6dSpjxoxh9erVfX7/CvQiEnl/+9vfuPLKK9mwYQNHHXUUDz/8cPLcm2++ydKlS9mwYQPNzc3ceOONTJkyhfPPP5+f/OQnNDU18elPf5pLLrmEW2+9lebmZioqKvjhD38IwGWXXcbixYt58cUXGThwYIef++KLL3L//ffzzDPPUFJSwtKlS3nppZd49tln+e53v0v7ygSbN2/m29/+Ns3Nzbz66qs0NDTw/PPPc9ttt3HLLbf0+f0r0ItI5I0ePZpYLAbA5z//ebZt25Y8d8QRR1BSUsLll1/OI488wuGHH35I+z179vD2228zbdo0AC699FJWrlzJ22+/TWtrK1OmTAGgtra2Q7sZM2bwiU98AoivYnnDDTdQWVnJmWeeyeuvv86//vWvZP8qKioYMGAAJ598MmeccQZmRkVFRYe+9laoQG9mNWa2ycw2m9n8NOdnmVmzmTWZ2Voz+0LYtiIi2faxj30s+XjgwIEcOHAg+XzQoEGsXr2ar3zlKzz66KPU1NSEft3u1gprX/4Y4IEHHmDXrl2sW7eOpqYmjj32WNra2g7p34ABA5LPBwwY0KGvvdVtoDezgcBdwFlAOXCRmaXenfgDcIq7x4BvAPf0oK2ISM7s3buXPXv2cPbZZ3P77bcnN0cJLjd85JFHcvTRR/Pcc88B8Ktf/Ypp06Zx9NFHU1ZWxqpVqwB48MEHO/05e/bsYfjw4QwePJhnn32Wf/zjH9l9YwFhbgNPBDa7+xYAM3sQmAVsbL/A3fcGri8FPGxbEZFcam1tZdasWbS1teHu/OxnPwNgzpw5XHHFFdx555089NBD3H///cybN4/33nuPE088kV/+8pcA/OIXv+CKK66gtLSU6upqjjzyyLQ/5+KLL+a8886jqqqKWCzG2LFj++09drtMsZn9B1Dj7pcnnn8NOM3dr0q57gLgx8Bw4Bx3fzFs21Rapjh6sjlyRmPq81chLVPcW3v37uXjH/84AAsXLuSNN97gjjvuyOrPzMYyxen2ujrk28Hdl7r7WODLwI960jbRybpEfX/trl27QnRLRCT3Hn/8cWKxGOPGjeO5557jxhtvzHWXDhGmdNMCHB94PgLY3tnF7r7SzD5tZsf0pK271wP1EM/oQ/RLRCTnZs+ezezZs3PdjS6FyejXAGPMbLSZDQHmAMuCF5jZZyyxy62ZjQeGALvDtBURkezqNqN39wNmdhXwJDAQuNfdN5jZvMT5xcBXgEvM7ANgHzDb48X/tG2z9F5EJA+5O4k8UDKgN9u/hlp8wd2XA8tTji0OPL4VuDVsW4kmrTkvqUpKSti9ezdDhw5VsM8Ad2f37t2UlJT0qJ0WNRORrBkxYgQtLS1ogEXmlJSUMGLEiB61UaAXkawZPHgwo0ePznU3ip7WuhERiTgFehGRiFPpRrJCa8iL5A9l9CIiEadALyIScSrdSM60Nj3BuxsbQ11bWl5NWSz8OuEicpAyesmZdzc2sn/n1m6v279za+gvBBE5lDJ66VfBLH7/zq0MGT6aT9Yu7LLNjob57N+5lR0NBzcoU4YvEp4CvfSr9ix+yPDRDBk+mtLy6m7bpF7T/ltAZ4FeSzGIdKRAL/0uTBYfVBar6RDUg5m9iHRPgV6yLl25pq+CpZz6E1+nrq6uz68pElUK9JJ1vSnXdCXYfv/OrVyzYBG3bDmuT6+pco9EmQK99Iuelmu6EizlqIwj0j0Feil4GpEj0jUFeumTdGvapE6EylRdPp2ejsgRKUYK9JJxwZo8kJG6fGc0Ikekewr0khWZrMmLSN8o0EvkaOilSEda60YipbS8Olky2r9zKw0NDTnukUjuKaOXjMjGpKje0NBLkUOFyujNrMbMNpnZZjM75P8eM7vYzJoTf14ws1MC57aZ2ctm1mRmazPZeckfwZUos3nzVUR6rtuM3swGAncBM4AWYI2ZLXP3jYHLtgLT3P0tMzsLqAdOC5yf7u7/zmC/JQ/pBqxIfgpTupkIbHb3LQBm9iAwC0gGend/IXD9KmBEJjsp+aPQ9oJtamqiuroagNraWt2YlaIUpnRzHPBa4HlL4lhn/hP4XeC5A0+Z2Toz0/9l0m9Ky6uJxWJAPODrxqwUqzAZvaU55mkvNJtOPNB/IXB4qrtvN7PhwO/N7FV3X5mmbR1QBzBy5MgQ3RLpWlmshsaFiwCSWb1IMQoT6FuA4wPPRwDbUy8ys0rgHuAsd9/dftzdtyf+3mlmS4mXgg4J9O5eT7y2T1VVVdovEskv+TLSRkS6FqZ0swYYY2ajzWwIMAdYFrzAzEYCjwBfc/e/Bo6XmllZ+2NgJrA+U52X3NJIG5HC0G1G7+4HzOwq4ElgIHCvu28ws3mJ84uBm4ChwH+bGcABd68CjgWWJo4NAhrc/YmsvBPJiUIaabNqy26tOy9FKdSEKXdfDixPObY48Phy4PI07bYAp6QeF8mF1OWMtTyCFAstgSCRNmr+44ya/zgvl1R0uIeg5RGkmGgJBOlWoY2dT0fLGUsxU0YvIhJxyugltP7cOUpEMkcZvYQWHE4JhT+kcuWqtZSMrKBkZAWtTRoMJtGljF56pJCGU3Yl+AWlfWYl6hTopShp3XopJirdiIhEnAK9iEjEKdCLiEScavSSVhQmSfVEcHkELY0gUaNAL10qhqWIU0fgNDQ0KNBLpCjQS5fax84PGT664MfNd0YjcCTqFOgF6LpUE5Wx8yLFSjdjRUQiThm9SDe0WYkUOmX0IiIRp4xeJEVTUxPV1dXJ560lFVoHRwqaAr1IQGl5Ne9ubGTVlt1A+5DS3Qr0UtAU6EUCtBOVRJECvRyiGCZJiRQT3YyVQwQ3GInqJCmRYhIqozezGuAOYCBwj7svTDl/MfC9xNO9wDfd/S9h2kp+0iQpkejoNtCb2UDgLmAG0AKsMbNl7r4xcNlWYJq7v2VmZwH1wGkh24rkNS14JoUuTOlmIrDZ3be4+37gQWBW8AJ3f8Hd30o8XQWMCNtWJJ+Vllcn71G0L3gmUmjClG6OA14LPG8BTuvi+v8EftfLtiJ5RQueSRSECfSW5pinvdBsOvFA/4VetK0D6gBGjhwZolsiIhJGmNJNC3B84PkIYHvqRWZWCdwDzHL33T1pC+Du9e5e5e5Vw4YNC9N3EREJIUygXwOMMbPRZjYEmAMsC15gZiOBR4Cvuftfe9JWRESyq9vSjbsfMLOrgCeJD5G81903mNm8xPnFwE3AUOC/zQzgQCI7T9s2S+9FJOuC6+DU1tZqBI4UhFDj6N19ObA85djiwOPLgcvDtpX8Ul9fz46GRcnnmg2bXml5NRVtLwPxgA8o0EtB0MxYoaGhITkTFjQbtjNlsRoaGxtpbGwkFovlujsioWmtGwE0Ezas9k1IdmzZzaQTh+a4NyLhKKMXEYk4BXoRkYhToBcRiTjV6EV6KXXLQQ23lHylQC/SC8GhlqDhlpLfFOhFeqEsVkPjwoNzD4KZvUi+UY1eRCTilNEXqfr6+uTa6k1NTXDE8V03EJGCpUBfpBoaGmhqaiIWixGLxXi5pCLXXSo47ZOnQBOoJL8p0BexWCxGY2Mj0DFoSe9owTPJVwr0IhmgBc8knynQi2RAcBSORuBIvtGoGxGRiFNGX8RWbdmt2nwGaWVLyVfK6EVEIk6BXiQL2kfgVFdXU19fn+vuSJFT6UYkwzQCR/KNAn0R0WzY/qEROJJvVLopIu2zYSE+WUr7wooUh1CB3sxqzGyTmW02s/lpzo81sxfN7H0zuy7l3DYze9nMmsxsbaY6Lr3TPhu2sbGRslhNrrsjIv2g29KNmQ0E7gJmAC3AGjNb5u4bA5e9CVwDfLmTl5nu7v/uY19FClLqMNZtC8/JYW+kGIWp0U8ENrv7FgAzexCYBSQDvbvvBHaamf4Fi6TYv3MrOxoO/iJcf+Lrujkr/SpM6eY44LXA85bEsbAceMrM1pmZ/nXnWHt2qYlS/aO2tpYhw0cnn+/fuTV5Q1ykv4TJ6C3NMe/Bz5jq7tvNbDjwezN71d1XHvJD4l8CdQAjR47swcuL5K+6ujpu2XIwLwpm9iL9JUxG3wIEx+GNALaH/QHuvj3x905gKfFSULrr6t29yt2rhg0bFvblRUSkG2EC/RpgjJmNNrMhwBxgWZgXN7NSMytrfwzMBNb3trMiItJz3ZZu3P2AmV0FPAkMBO519w1mNi9xfrGZfRJYCxwBfGRm1wLlwDHAUjNr/1kN7v5EVt6JiIikFWpmrLsvB5anHFsceLyDeEkn1TvAKX3poPRecCYsaDZsvtBOVNLfNDM2woIzYUGzYfNBaXk1sVgMiAd8jcCR/qC1biIuuC8saG/YXNM6OJILCvQRp81FRESlGxGRiFNGL5IFYX+L0o1Z6Q8K9CI5Ultbm3ysDUokmxToRXKkrq4uGdh1Y1aySTV6EZGIU6AXEYk4lW4iRvvCFq7gjVnQzVnJHGX0EaN9YQtTbW1tcsYsaNasZJYy+ggKzobVZKnCkLpufdsWrVsvmaNAL5JDXX0RB2c1a59Z6QuVbkREIk4ZvUg/C1tOC24qrg3FpS+U0YvkodLy6uSm4tpQXPpKGb1IHiqL1VAWqwG0obj0nTJ6EZGIU6CPgPr6eqqrq6muru6wo5SICCjQR0LqJKngqogSDStXraVkZAUlIyuor6/PdXekwKhGHxGpWwZKdARnN7ffmNUIHOkJBXqRPKcbs9JXoUo3ZlZjZpvMbLOZHfIvzczGmtmLZva+mV3Xk7aSeaPmP578IyLSbUZvZgOBu4AZQAuwxsyWufvGwGVvAtcAX+5FWxHpAa1yKT0VJqOfCGx29y3uvh94EJgVvMDdd7r7GuCDnrYVkfBKy6u1yqX0WJga/XHAa4HnLcBpIV+/L21FJEVZrIbGhYuSz7UFoYQRJqO3NMc85OuHbmtmdWa21szW7tq1K+TLi4hId8IE+hYguE3RCGB7yNcP3dbd6929yt2rhg0bFvLlRUSkO2FKN2uAMWY2GngdmAOEnZHTl7bSieB2gRCv0wbrtlJcgjdndWNW0uk20Lv7ATO7CngSGAjc6+4bzGxe4vxiM/sksBY4AvjIzK4Fyt39nXRts/Reikb7TNj24K7ZsMUr+N+9fXa0Ar2kCjVhyt2XA8tTji0OPN5BvCwTqq30nWbCCsSDentg141Z6YxmxopEVOqEOW1HWLwU6EUiJFivby2pSC6dIMVNgV4kIlLr9W1H7FagF0DLFItERl1dHY2NjTQ2NmoUlnSgjF4kooKbi4M2GC9myugLhHaRkp6ora1Nbi4O2mC82CmjLxDBsfMaN1/culp+un1kTV1dHbdsOS55XOvYFzcF+gKisfMi0hsq3YgUifahl9XV1dp3tsgooxcpAqXl1VS0vQxoqYRipEAvUgTKYjVsIz6mvm2L6vXFRoFeJELC7hOs7QiLiwJ9RGgjcAkrWMYBlXKKgQJ9HguuO6815yVTgmUcUCmnGGjUTR5rHzsPWnNesksjcqJNGX2e09h5yTaNyIk+BXqRIlcWq6Fx4SJAm5dElQJ9gdLNV8kW7UEbPQr0IpKkPWijSYG+gCiLl2xL3YNW4+2jQYE+jwSHU4KGVEpupY7yUoZfuBTo80hwKWLQkErJrWB2D7pRW8hCBXozqwHuAAYC97j7wpTzljh/NvAe8HV3fylxbhvQCnwIHHD3qoz1PoI0nFJyIUxZcMeW3ezfuZWSkRUA3Hnj1cruC0S3E6bMbCBwF3AWUA5cZGblKZedBYxJ/KkD/ifl/HR3jynIixSu0vLq5K5V2rGqsITJ6CcCm919C4CZPQjMAjYGrpkF/K+7O7DKzI4ys0+5+xsZ77GI5ERZrIayWHzphB0N8zUMs4CEWQLhOOC1wPOWxLGw1zjwlJmtMzP9SxCJgNLy6uS9pKamJmX3eS5MRm9pjnkPrpnq7tvNbDjwezN71d1XHvJD4l8CdQAjR44M0S0RyZXU2bQahpnfwmT0LcDxgecjgO1hr3H39r93AkuJl4IO4e717l7l7lXDhg0L1/sIqK+vTy4m1T58TaSQ1NbWdhgGrAw//4QJ9GuAMWY22syGAHOAZSnXLAMusbhJwB53f8PMSs2sDMDMSoGZwPoM9r/gaYVKKXR1dXU0NjYm/8RiMa2GmWe6Ld24+wEzuwp4kvjwynvdfYOZzUucXwwsJz60cjPx4ZWXJZofCyyNj75kENDg7k9k/F0UOA2plEIUHJK5beE5ycdaRiH/hBpH7+7LiQfz4LHFgccOXJmm3RbglD72UUQKSFfLKKh2nxuaGSsiWRPM7lesWMGKFSs61O8V+PuHAn0/03o2UgwOlnWOY1uiLJnu3z6orNMfFOj7mdazkSgKs4RCXV0dt2w5OAWnbYsmXfUXBfoc0M1XkY5bGKaWdRT0M0uBPs9pDXqJquCkq2BZRyWdzFOgF5GcCdbymXQ9cGhJB5Th95UCfT9IzVZ081WKSU9/Kw2WdOBgWeeaBfHsX8sj91yYmbHSR5r9KhJeWaymw0zbu+++m48dPw6A919bz9y5czXrtocsPtcpv1RVVfnatWtz3Y2Maf8VtDc3YFWjFzmotemJDjdwAaZNm5Y8X8wlHjNb19meH8roRaRgBLP9u+++u0OQX7FiBXPnzqVkZAUlIyuU7QeoRp8lva3LK4MXCSd1XP4nSip4d2MjcLDEo+GacQr0WRKcGKW6vEj2BXfASi3xFPvSCwr0WRR2YpSyeJHMSh2jf82CRazashuIZ/upgT8oil8CCvQZ0pM1bBTYRXqvp///pJZ4gtl+qq5m6Kb+3ODSzPlOo24ypH051mBw7ywzUKAXya3OgnR79g/xzB9IDu1MVVpenSwVhXnt7nS2vn9YXY26UUafQVrDRqSwBbP/1qYnkjd3U73/2nref2192vP1J77eaemnr8G8txTo+0AzXkWiK3hzN1VnXwKpo30OaVdS0elrZpMCfR9oZI1IYepr+bSzL4Hu6v+wIvkFUb3qJx3O70jcLIaufyvoDdXoeyhdFq+RNSLSLkz9vyvt9waCk8E6E7wPqBp9BimLF5HeSB3905nU3wpWBTJ9gEknDgXSzw/ojAJ9NzobNqmbriKSTibKQu1zANK93rbE358oqei0TJRKgT6NYHBPXTipqyy+kMfZikj+CPNlURarYRuB+wSJWJVOqEBvZjXAHcBA4B53X5hy3hLnzwbeA77u7i+FaZsvOgvu06ZNi+RMOREpHt0GejMbCNwFzABagDVmtszdNwYuOwsYk/hzGvA/wGkh2+ZEakmmt8G9q29e3YAVkXwQJqOfCGx29y0AZvYgMAsIButZwP96fAjPKjM7ysw+BYwK0bbPUoN2GKklGWXuIhJVYQL9ccBrgectxLP27q45LmTbQ2zatKnDfpHdSbcBQXf6krWr9i4ihSRMoLc0x1IH33d2TZi28RcwqwPao+7eFStWbArRtw5WdHEzorPr586d29npY4B/pztht/asXxHQ6WdRhPRZHKTP4qB8+CxO6OxEmEDfAhwfeD4C2B7ymiEh2gLg7vVA3mwJY2ZrO5t8UGz0WRykz+IgfRYH5ftnEWYrwTXAGDMbbWZDgDnAspRrlgGXWNwkYI+7vxGyrYiIZFG3Gb27HzCzq4AniQ+RvNfdN5jZvMT5xcBy4kMrNxMfXnlZV22z8k5ERCStUOPo3X058WAePLY48NiBK8O2LRB5U0bKA/osDtJncZA+i4Py+rPIy0XNREQkc8LU6EVEpIAp0IdgZteZmZvZMbnuS66Y2U/M7FUzazazpWZ2VK771J/MrMbMNpnZZjObn+v+5IqZHW9mz5rZK2a2wcy+nes+5ZqZDTSzP5vZY7nuS2cU6LthZscTX8Lhn7nuS479Hhjn7pXAX4Hv57g//SawlMdZQDlwkZmV57ZXOXMA+K67fw6YBFxZxJ9Fu28Dr+S6E11RoO/ez4D/RycTvYqFuz/l7gcST1cRnxNRLJLLgLj7fqB9KY+i4+5vtC9Y6O6txANc94usR5SZjQDOAe7JdV+6okDfBTM7H3jd3f+S677kmW8Av8t1J/pRZ0t8FDUzGwWcCvwpx13JpduJJ4If5bgfXSr69ejN7Gngk2lO/RdwAzCzf3uUO119Fu7+28Q1/0X81/cH+rNvORZ6KY9iYWYfBx4GrnX3d3Ldn1wws3OBne6+zsyqc9ydLhV9oHf3M9MdN7MKYDTwl/hy+4wAXjKzie6+ox+72G86+yzamdmlwLnAGV5c43LDLANSNMxsMPEg/4C7P5Lr/uTQVOB8MzsbKAGOMLNfu/tXc9yvQ2gcfUhmtg2ocvdcL1yUE4kNZH4KTHP3XbnuT38ys0HEb0CfAbxOfGmP2mKc5Z3YZOh+4E13vzbH3ckbiYz+Onc/N8ddSUs1egnr50AZ8HszazKzxd01iIrETej2pTxeAf6vGIN8wlTga8AXE/8OmhIZreQxZfQiIhGnjF5EJOIU6EVEIk6BXkQk4hToRUQiToFeRCTiFOhFRCJOgV5EJOIU6EVEIu7/A6RLZJZRm20nAAAAAElFTkSuQmCC\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}(x)f_Y(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}(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}(x)f_Y(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}\n",
|
|
"$$\n",
|
|
"using that $r^2 = x^2 + y^2$ in the last step. We conclude that $X$ and $Y$ are independent, and 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": "code",
|
|
"execution_count": null,
|
|
"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": [],
|
|
"source": [
|
|
"def random_normal_pair():\n",
|
|
" '''Return two independent normal random variables.'''\n",
|
|
" # YOUR CODE HERE\n",
|
|
" raise NotImplementedError()\n",
|
|
" return x, y\n",
|
|
"\n",
|
|
"# Plotting\n",
|
|
"# YOUR CODE HERE\n",
|
|
"raise NotImplementedError()"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python 3",
|
|
"language": "python",
|
|
"name": "python3"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": {
|
|
"name": "ipython",
|
|
"version": 3
|
|
},
|
|
"file_extension": ".py",
|
|
"mimetype": "text/x-python",
|
|
"name": "python",
|
|
"nbconvert_exporter": "python",
|
|
"pygments_lexer": "ipython3",
|
|
"version": "3.9.12"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 5
|
|
}
|