mp573_hw4_handout.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework Set 4 - handout\n",
    "\n",
    "### Please see computational questions below, and comments labeled 'To-do', which indicate the code that needs to be inserted by you"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "## MP 573 HW4\n",
    "##\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import display, Image\n",
    "import matplotlib.image as mpimg\n",
    "from os.path import dirname, join as pjoin\n",
    "import scipy.io as sio\n",
    "from scipy import signal, ndimage, misc\n",
    "from skimage.data import shepp_logan_phantom\n",
    "import scipy.stats as stats\n",
    "\n",
    "font = {'weight' : 'normal',\n",
    "        'size'   : 20}\n",
    "\n",
    "np.set_printoptions(precision=2)\n",
    "np.set_printoptions(suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 1.1: Gaussian vs  Bernoulli \n",
    "\n",
    "Suppose we have two independent random variables ($X$ and $Y$),  one of which ($X$) is Gaussian distributed with mean $\\mu=0$ and standard deviation $\\sigma=1$,  and the other ($Y$) is Bernoulli distributed with $p=0.5$.  What is the probability that $X>Y$? Please report your empirical and your theoretical results,  and show your work/logic for your theoretical results.  Do the empirical results agree with your theoretical calculation? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Empirical probability Gaussian > Bernoulli = 0.32815\n"
     ]
    }
   ],
   "source": [
    "# Number of trials\n",
    "N = 100000;\n",
    "\n",
    "# Distribution parameters (Bernoulli)\n",
    "p = 0.5;\n",
    "\n",
    "# Distribution parameters (Gaussian)\n",
    "mu = 0\n",
    "sigma = 1\n",
    "\n",
    "# Generate some (pseudo-)random data\n",
    "y = stats.bernoulli.rvs(p, size=N) # Bernoulli\n",
    "x = stats.norm.rvs(mu, sigma, size=N)\n",
    "\n",
    "# Count the number of instances where the Gaussian data were larger than the Bernoulli data\n",
    "Ne = sum(x>y)\n",
    "\n",
    "# Calculate the empirical fraction of experiments where the Gaussian (mean=0, std=1) was larger than the Bernoulli (p=0.5)\n",
    "EmpProb = Ne/N\n",
    "print('Empirical probability Gaussian > Bernoulli = ' + str(EmpProb))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 1.2: Magnitude of Gaussian vs Bernoulli \n",
    "\n",
    "As before,  suppose we have two independent random variables ($X$ and $Y$),  one of which ($X$) is Gaussian distributed with mean $\\mu=0$ and standard deviation $\\sigma=1$,  and the other ($Y$) is Bernoulli distributed with $p=0.5$.  Now suppose we generate a third random variable $Z = |X|$ (i.e.,  $Z$ is the magnitude of our Gaussian random variable).   \\emph{What is the probability that $Z>Y$?} Please report your empirical and your theoretical results,  and show your work/logic for your theoretical results.  Do the empirical results agree with your theoretical calculation? \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Empirical probability magnitude of Gaussian > Bernoulli = 0.6541\n"
     ]
    }
   ],
   "source": [
    "# Number of trials\n",
    "N = 10000;\n",
    "\n",
    "# Distribution parameters (Bernoulli)\n",
    "p = 0.5;\n",
    "\n",
    "# Distribution parameters (Gaussian)\n",
    "mu = 0\n",
    "sigma = 1\n",
    "\n",
    "# Generate some (pseudo-)random data\n",
    "y = stats.bernoulli.rvs(p, size=N) # Bernoulli\n",
    "x = stats.norm.rvs(mu, sigma, size=N)\n",
    "z = abs(x)\n",
    "\n",
    "# Count the number of instances where the Gaussian data were larger than the Bernoulli data\n",
    "Ne = sum(z>y)\n",
    "\n",
    "# Calculate the empirical fraction of experiments where the Gaussian (mean=0, std=1) was larger than the Bernoulli (p=0.5)\n",
    "EmpProb = Ne/N\n",
    "print('Empirical probability magnitude of Gaussian > Bernoulli = ' + str(EmpProb))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 1.3: Effect of the number of trials\n",
    "\n",
    "For the previous question (Magnitude of Gaussian vs Bernoulli), measure $N_e/N$ empirically as described above, using varying numbers of trials $N$. Specifically, use the following values of $N$: $20$, $10^3$, $10^5$. For each of these values of $N$, repeat the experiment at least 10 times, leading to at least 10 empirical measurements of $N_e/N$. For example, for $N=20$, you would run the code to generate random numbers following the distributions of $Y$ and $Z$ (20 each), leading to an estimate of $N_e/N$, and would repeat this process at least 10 times (leading to at least 10 different estimates of $N_e/N$) before moving on to the next value of $N$. For each value of $N$, report (or plot as a function of $N$), the empirical average of $N_e/N$ for that value of $N$, as well as the empirical standard deviation of $N_e/N$ for that value of $N$. If you plot the results, use  $\\log(N)$ as you x axis, so that the results can be observed more easily. \n",
    "\n",
    "Describe (with words)  your observations of the behavior of the mean and standard deviation of $N_e/N$ as you increase $N$. \n",
    "\n",
    "Feel free to complete this problem manually by clicking run multiple times on the corresponding Jupyter notebook cell and writing down the results, or automatically by writing a loop that will automatically run this code multiple times and record the numbers for you. The second approach is recommended for reproducibility and minimization of manual work, but either approach is valid for this homework set. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To-do: \n",
    "# For a given value of N, obtain multiple (at least 10) empirical measurements \n",
    "# of Ne/N based on the code above\n",
    "# And repeat this process for multiple values of N, in order to estimate and plot the \n",
    "# mean and standard deviation of Ne/N as a function of N. \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 2.3: Sum of independent random variables \n",
    "\n",
    "### General case \n",
    "\n",
    "Suppose we have two independent random variables, $X$ and $Y$, with probability density functions $f_X(x)$ and $f_Y(x)$, respectively. Derive the probability density function of $f_Z(x)$ of $Z = X+Y$, in terms of $f_X(x)$ and $f_Y(x)$. \n",
    "\n",
    "\n",
    "### Sum of two uniform random variables\n",
    "\n",
    "Suppose that both $X$ and $Y$ are uniformly distributed random variables between $a=-0.5$ and $b=0.5$. What is the distribution of $Z = X+Y$? \n",
    "\n",
    "### Computational confirmation (10 points)\n",
    "\n",
    "Confirm your result from the previous subsection (sum of two uniform random variables) computationally. You may consider  using code similar to that provided at the end of the Jupyter notebook for lecture 20. If you want to show a histogram normalized as a probability density function (i.e., so it integrates to 1), you may consider using the hist function (as in lecture 20), adding the option 'density=True'. \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To-do: Insert code here to generate data sampled from a random variable\n",
    "# that is the sum of two independent, uniformly-distrbuted random variables. \n",
    "# You may want to check the code from lecture 20 for this purpose\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 3.2: Linear transformation of Gaussian noise\n",
    "\n",
    "### Theory \n",
    "\n",
    "Suppose we have a vector random variable $\\mathbf{X}$ consisting of $M$ components $X_m$, each of which is an independent Gaussian random variable with mean $\\mu=0$ and variance $\\sigma^2=1$. Next, we generate a new random variable $\\mathbf{Y}$ by linear transformation of $\\mathbf{X}$, as $\\mathbf{Y} = \\mathbf{A X}$ where $\\mathbf{A}$ is some $M \\times M$ real-valued matrix. Note that this example is representative of noise propagation in imaging, as the noise propagates from the measurements to the reconstructed image. \n",
    "\n",
    "What is the covariance matrix of $\\mathbf{X}$?\n",
    "\n",
    "What are the mean vector and covariance matrix of $\\mathbf{Y}$?\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### Computational example\n",
    "\n",
    "Now consider the specific case when $M=3$ (i.e., $\\mathbf{X}$ has three i.i.d. components), and \n",
    "$$\\mathbf{A} = \\left(  \\begin{array}{ccc}   1& 2& 3 \\\\ 2 & 1& 2 \\\\1 & 0 & 1 \\end{array}    \\right)$$\n",
    "\n",
    "Calculate a large number (e.g., $N=10^5$) of vectors $\\mathbf{x}$, distributed according to the probability distribution of $\\mathbf{X}$, and generate the corresponding vectors $\\mathbf{y} = \\mathbf{A x}$. Calculate the empirical covariance matrix for $\\mathbf{X}$ and $\\mathbf{Y}$, and also calculate the theoretical covariance matrices for $\\mathbf{X}$ and $\\mathbf{Y}$. Compare the theoretical and empirical results. Please see the associated Jupyter notebook for some code to get started. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " Covariance matrix of X (empirical):\n",
      "[[ 1.    0.01  0.01]\n",
      " [ 0.01  0.99 -0.03]\n",
      " [ 0.01 -0.03  1.02]]\n",
      "\n",
      " Covariance matrix of X (theoretical):\n",
      "[[0. 0. 0.]\n",
      " [0. 0. 0.]\n",
      " [0. 0. 0.]]\n",
      "\n",
      " Covariance matrix of Y (empirical):\n",
      "[[13.93 10.03  4.06]\n",
      " [10.03  9.07  4.06]\n",
      " [ 4.06  4.06  2.04]]\n",
      "\n",
      " Covariance matrix of Y (theoretical):\n",
      "[[0. 0. 0.]\n",
      " [0. 0. 0.]\n",
      " [0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "# Generate a bunch of Gaussian iid 3x1 vectors (X)\n",
    "N = 10000\n",
    "mu = 0\n",
    "sigma = 1\n",
    "x = stats.norm.rvs(mu, sigma, size=[3,N])\n",
    "\n",
    "# Generate some matrix A\n",
    "A =  np.array([[1, 2, 3], [2, 1, 2], [1, 0, 1]])\n",
    "\n",
    "# Generate our transformed random data Y=AX\n",
    "y = np.matmul(A,x)\n",
    "\n",
    "# Calculate means\n",
    "mux = np.mean(x,axis=1)\n",
    "muy = np.mean(y,axis=1)\n",
    "\n",
    "# Calculate empirical covariance matrices for X and Y\n",
    "Cx = np.cov(x)\n",
    "Cy = np.cov(y)\n",
    "\n",
    "# Theoretical covariance matrices for X and Y\n",
    "# To-do by student: \n",
    "# Please replace the right-hand side of the following two lines\n",
    "# with the correct expression for the covariance matrix \n",
    "Cx_theory = np.zeros((3,3))\n",
    "Cy_theory = np.zeros((3,3))\n",
    "\n",
    "\n",
    "print('\\n Covariance matrix of X (empirical):')\n",
    "print(str(Cx))\n",
    "\n",
    "print('\\n Covariance matrix of X (theoretical, needs to be provided by student):')\n",
    "print(str(Cx_theory))\n",
    "\n",
    "print('\\n Covariance matrix of Y (empirical):')\n",
    "print(str(Cy))\n",
    "\n",
    "print('\\n Covariance matrix of Y (theoretical, needs to be provided by student):')\n",
    "print(str(Cy_theory))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}