homeworks/hw_02.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework 2\n",
    "\n",
    "**Due: 01/30/2018**\n",
    "\n",
    "## Instructions\n",
    "\n",
    "+ In any case, develop the code and generate the figures you need to solve the problems using this notebook.\n",
    "+ For the answers that require a mathematical proof or derivation you can either:\n",
    "    \n",
    "    - Type the answer using the built-in latex capabilities. In this case, simply export the notebook as a pdf and upload it on gradescope; or\n",
    "    - you can print the notebook (after you are done with all the code), write your answers by hand, scan, turn your response to a single pdf, and upload on gradescope.\n",
    "\n",
    "+ The total homework points are 100. Please note that the problems are not weighed equally."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Student details\n",
    "\n",
    "+ **First Name:**\n",
    "+ **Last Name:**\n",
    "+ **Email:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Here are some modules that you may need - please run this block of code:\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "sns.set_context('paper')\n",
    "import numpy as np\n",
    "import scipy.stats as st"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 1\n",
    "### Part A\n",
    "Let $X$ be a continuous random variable with CDF:\n",
    "$$\n",
    "F(x) = p(X \\le x).\n",
    "$$\n",
    "Show that the random variable\n",
    "$$\n",
    "Z = F(X)\n",
    "$$\n",
    "    is distributed uniformly in $[0,1]$.\n",
    "    Hint: Show that:\n",
    "$$\n",
    "    F_Z(z) := p(Z \\le z) = z.\n",
    "$$\n",
    "**Proof:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Type your proof here. Delete that ``<br>`` line (it just makes some white space).*\n",
    "<br><br><br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The theorem you have just proved is very useful when you want to test if a set of observations $x_1,x_2,\\dots,x_N$ has been indeed independent realizations of a random variable $X$ with CDF $F(x)$.\n",
    "Part A proves that, if this hypothesis is valid, then the transformed dataset $z_1,z_2,\\dots,z_N$, where\n",
    "$$\n",
    "z_i = F(x_i),\n",
    "$$\n",
    "should be distributed uniformely in $[0,1]$.\n",
    "In other words, the empirical histrogram of $z_1,z_2,\\dots,z_N$ should be a flat line.\n",
    "Use this observation to find the distribution from which this dataset was sampled:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = np.loadtxt('hw2_p1b_data.txt')\n",
    "fig, ax = plt.subplots()\n",
    "ax.hist(data, normed=True, alpha=0.5)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('Empirical PDF')\n",
    "ax.set_title('Histogram of data')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The correct distribution is one of the following:\n",
    "\n",
    "1. Standard normal, $\\mathcal{N}(0,1)$;\n",
    "\n",
    "2. Normal with mean 2 and variance 2, $\\mathcal{N}(2,2); or\n",
    "\n",
    "3. Exponential with rate parameter $1$, $\\mathcal{E}(1)$; or\n",
    "\n",
    "4. Exponential with rate parameter $2$, $\\mathcal{E}(2)$; or\n",
    "\n",
    "5. Exponential with rate parameter $10$, $\\mathcal{E}(10)$; or\n",
    "\n",
    "6. Gamma distribution with parameters $\\alpha=2.$ and $\\beta=3.$.\n",
    "\n",
    "Systematically go over these distributions and try to determine which one generated the data.\n",
    "All the required CDF's and inverse CDF's are implemented in [scipy.stats](https://docs.scipy.org/doc/scipy/reference/stats.html).\n",
    "Check also [scipy.stats.rv_continuous](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.html#scipy.stats.rv_continuous).)\n",
    "Please pay special attention to the defintiion of the probability distributions of the various random variables and how you can control their parameters.\n",
    "As a hint, here is how you can test for $\\mathcal{N}(0,1)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Testing for N(0,1):\n",
    "transformed_data = st.norm(loc=0, scale=1).cdf(data) # cdf() gives the CDF of a random variable\n",
    "# If the data came from N(0,1), the histogram of the transformed_data should match that of a uniform:\n",
    "fig, ax = plt.subplots()\n",
    "ax.hist(transformed_data, normed=True, alpha=0.5)\n",
    "ax.plot(np.linspace(0, 1,50), np.ones(50), lw=2, color='r') # This is the line that you should try to much.\n",
    "ax.set_xlabel('$z = F(x)$')\n",
    "ax.set_ylabel('Empirical PDF')\n",
    "ax.set_title('Testing the $\\mathcal{N}(0,1)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 2\n",
    "\n",
    "### Part A\n",
    "Let $U\\sim \\mathcal{U}(0,1)$ and let $F(x)$ be a function with the following properties:\n",
    "+ $F(-\\infty) = 0$,\n",
    "+ $F(\\infty)=1$, \n",
    "+ and $x_1 \\le x_2\\Rightarrow F(x_1)\\le F(x_2)$.\n",
    "\n",
    "Consider the random variable:\n",
    "$$\n",
    "X = F^{-1}(U),\n",
    "$$\n",
    "where $F^{-1}$ is the inverse of $F$.\n",
    "Show that the CDF of $X$ is $F(x)$, i.e.,\n",
    "$$\n",
    "F(x) = P(X \\le x).\n",
    "$$\n",
    "This is called *inverse transform sampling*.\n",
    "In this part (Part A), you have to prove it.\n",
    "In a while, you will use it to learn how to sample from a continuous univariate random variable when all you have is samples from it.\n",
    "\n",
    "**Proof:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Type your proof here. Delete that ``<br>`` line (it just makes some white space).*\n",
    "<br><br><br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part B\n",
    "\n",
    "Consider the following data set $x_1,\\dots,x_N$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = np.loadtxt('hw2_p2_data.txt')\n",
    "data = np.array(data)\n",
    "fig, ax = plt.subplots()\n",
    "ax.hist(data, alpha=0.5)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('Empirical PDF')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your goal is to generate a procedure that samples from the same distribution as the observed data. \n",
    "This is a variation of the standard *density estimation* problem.\n",
    "In general, this is a very difficult problem and we will see various ways to solve it later on.\n",
    "In this problem, you will develop a simple method that relies on the empirical CDF of the observed data.\n",
    "Needless to say, this method works only for one dimensional cases in which you have a lot of observations.\n",
    "\n",
    "The [empirical CDF](https://en.wikipedia.org/wiki/Empirical_distribution_function) of our data set, $x_1,\\dots,x_N$ is defined to be:\n",
    "$$\n",
    "\\hat{F}_N(x) = \\frac{\\text{Number of observations}\\;\\le x}{N} = \\frac{1}{N}\\sum_{i=1}^N 1_{[x_i,+\\infty]}(x),\n",
    "$$\n",
    "where $1_A(x)$ is the [indicator function](https://en.wikipedia.org/wiki/Indicator_function) of the set $A$.\n",
    "Using the, so called, [strong law of large numbers](https://en.wikipedia.org/wiki/Law_of_large_numbers#Strong_law), we can show that $\\hat{F}_N(x)$ converges to the true CDF of the data.\n",
    "\n",
    "Complete the code that calculates the empirical CDF:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def myECDF_base(x):\n",
    "    \"\"\"\n",
    "    Make this code work if ``x`` is a simple scalar.\n",
    "    \n",
    "    :param x:    The point at which you want to observe the PDF.\n",
    "    :param data: All the data you have observed (1D numpy array).\n",
    "    :returns:    The value of the empirical CDF at ``x``.\n",
    "    \"\"\"\n",
    "    N = data.shape[0]\n",
    "    # Write your code here (delete the next line and return the right value)\n",
    "    raise NotImplementedError('Implement me!')\n",
    "\n",
    "# Vectorize your function (i.e., make it work with 1D numpy arrays).\n",
    "# See this: https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.vectorize.html\n",
    "myECDF = np.vectorize(myECDF_base)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# You can test your results by comparing the empirical CDF you can get with of scipy.stats\n",
    "# The two should match almost exactly\n",
    "hist_rv = st.rv_histogram(np.histogram(data, bins=1000))\n",
    "fig, ax = plt.subplots()\n",
    "# The range in which the x's takes values:\n",
    "x_min = data.min()\n",
    "x_max = data.max()\n",
    "xx = np.linspace(x_min, x_max, 100)\n",
    "ax.plot(xx, myECDF(xx), label='My empirical CDF')\n",
    "ax.plot(xx, hist_rv.cdf(xx), '--', label='Empirical CDF from scipy.stats')\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('Empirical CDF')\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now complete the code that computes the inverse of the empirical CDF $\\hat{F}^{-1}$.\n",
    "There are may ways of doing this.\n",
    "Let's do it in a way that will teach us something about the root finding toolbox of numpy (see [this]()).\n",
    "Mathematically, we wish to find a function $F^{-1}$ such that\n",
    "$$\n",
    "F(F^{-1}(u))) = u,\n",
    "$$\n",
    "for any $u\\in[0,1]$ (the domain in which $F(x)$ takes values).\n",
    "It is obvious that $F^{-1}(u)$ is the solution to the *root finding* problem:\n",
    "$$\n",
    "F(x^*) = u.\n",
    "$$\n",
    "Since we know that $F$ is increasing, this problem must have a unique solution for any $u\\in[0,1]$.\n",
    "To find this solution, we can use [Brent's method](https://en.wikipedia.org/wiki/Brent%27s_method).\n",
    "Please note, that the problem that this code solves is of the form:\n",
    "$$\n",
    "g(x^*) = 0.\n",
    "$$\n",
    "So, you will have to reformulate the original problem as:\n",
    "$$\n",
    "F(x^*) - u = 0.\n",
    "$$\n",
    "Study the [numpy implementation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brentq.html#scipy.optimize.brentq) of Brent's method and complete the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy import optimize    # Gives you access to optimize.brentq\n",
    "\n",
    "\n",
    "\n",
    "def myiECDF_base(u):\n",
    "    \"\"\"\n",
    "    Evaluates the inverse of the empirical CDF.\n",
    "    \n",
    "    :param u:   A scalar at which to evaluate the function.\n",
    "    :param data: All the data you have observed (1D numpy array).\n",
    "    :returns:    The value of the inverse of the empirical CDF at ``u``.\n",
    "    \"\"\"\n",
    "    # Write your code here\n",
    "    # You will have to define a functioin to pass to optimize.brentq\n",
    "    # You can define this function in here\n",
    "    # (delete the next line and return the right value)\n",
    "    raise NotImplementedError('Implement me!')\n",
    "\n",
    "# Vectorize your function (i.e., make it work with 1D numpy arrays).\n",
    "# See this: https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.vectorize.html\n",
    "myiECDF = np.vectorize(myiECDF_base)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# You can test your results by comparing the inverse of the empirical CDF you can get with scipy.stats\n",
    "# The two should match almost exactly\n",
    "hist_rv = st.rv_histogram(np.histogram(data, bins=1000))\n",
    "fig, ax = plt.subplots()\n",
    "uu = np.linspace(1e-3, 1, 100)    # For convergence issues we cannot start at 0\n",
    "ax.plot(uu, myiECDF(uu), label='Inverse of my empirical CDF')\n",
    "ax.plot(uu, hist_rv.ppf(uu), '--', label='Inverse of empirical CDF from scipy.stats')\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('Inverse of empirical CDF')\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part C\n",
    "\n",
    "Now use the *inverse transform sampling* method to generate samples from same distribution as the original data.\n",
    "That is, you can now generate uniform samples:\n",
    "$$\n",
    "u_i \\sim U([0,1]),\n",
    "$$\n",
    "and transform them as:\n",
    "$$\n",
    "\\hat{x}_i = {\\hat{F}}^{-1}(u_i).\n",
    "$$\n",
    "The $\\hat{x}_i$'s generated in this way should have the same distribution of the data you started with.\n",
    "Verify this by comparing the histrogram of $1,000$ $\\hat{x}_i$ samples with the original data of this problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Write your code here.\n",
    "# Feel free to copy paste code from above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 3\n",
    "Demonstate that the central limit theorem holds for the data-generated distribution of Problem 2.\n",
    "Feel free to copy paste code from [handout 5](http://localhost:8888/notebooks/uq-course/handouts/handout_05.ipynb).\n",
    "At the very least you should\n",
    "\n",
    "1. Consider the sum:\n",
    "$$\n",
    "S_M = \\frac{X_1+\\dots+X_M}{M} = \\frac{\\hat{F}^{-1}(U_1)+\\dots+\\hat{F}^{-1}{U_M}}{M},\n",
    "$$\n",
    "and plot its histogram for various M;\n",
    "\n",
    "2. You should calculate $\\mu$ and $\\frac{\\sigma^2}{M}$ from the right-hand-side of the central limit theorem and draw the corresponding Gaussian along with the histogram for 1;\n",
    "\n",
    "3. You should repeat 1 and 2 for various values of $M$ until you start getting convergent results.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Put your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 4\n",
    "\n",
    "### Part A\n",
    "Let $X$ be a random variable and assume that you know that\n",
    "\n",
    "+ it's taking values in $\\mathbb{R}$,\n",
    "\n",
    "+ its mean is:\n",
    "$$\n",
    "\\mathbb{E}[X] = \\mu,\n",
    "$$\n",
    "\n",
    "+ and its variance is:\n",
    "$$\n",
    "\\mathbb{V}[X] = \\sigma^2.\n",
    "$$\n",
    "\n",
    "Using the principle of maximum entropy, show that\n",
    "$$\n",
    "X \\sim \\mathcal{N}(\\mu, \\sigma^2).\n",
    "$$\n",
    "\n",
    "**Proof:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Type your proof here. Delete that ``<br>`` line (it just makes some white space).*\n",
    "<br><br><br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part B\n",
    "Let $X$ be a random variable and assume that you know that\n",
    "\n",
    "+ it's taking values in $[0,\\infty)$, and\n",
    "+ it's mean is:\n",
    "$$\n",
    "\\mathbb{E}[X] = \\frac{1}{\\lambda}.\n",
    "$$\n",
    "\n",
    "Using the principle of maximum entropy, show that\n",
    "$$\n",
    "X\\sim\\mathcal{E}(\\lambda).\n",
    "$$\n",
    "\n",
    "**Proof:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Type your proof here. Delete that ``<br>`` line (it just makes some white space).*\n",
    "<br><br><br><br><br><br><br><br><br><br>"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}