{ "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": [ "
" ] }, "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": [ "
" ] }, "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": null, "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": [], "source": [ "def f_inv_discrete(prob_list,p):\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", "\n", "# plotting\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "code", "execution_count": null, "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": null, "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": [], "source": [ "def sample_Zn(alpha,b,n):\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", "\n", "# Plotting\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "code", "execution_count": null, "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": [ "YOUR ANSWER HERE" ] }, { "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": null, "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": [], "source": [ "from scipy.stats import levy_stable\n", "\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "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": [ "YOUR ANSWER HERE" ] }, { "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 }