Homework-4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework #4:  Nicolas Acton, ECE-5714"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the Gaussian mixture defined as:\n",
    "$$ G = (1 - \\varepsilon) N(0,1) + \\varepsilon N(0,25) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1\n",
    "Express the asymptotic variance of the normalized sample mean at G as a function of $\\varepsilon$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ f(r) = (1 - \\varepsilon)\\phi(r) + \\varepsilon\\phi(\\frac{r}{a}) = \\frac{(1 - \\varepsilon)}{\\sqrt{2\\pi}}e^{-\\frac{r^{2}}{2}} + \\frac{(\\varepsilon)}{\\sqrt{2\\pi}a}e^{-\\frac{r^{2}}{2a^{2}}}$$\n",
    "\n",
    "$$ lim_{m\\rightarrow\\infty}V(\\sqrt{m}\\bar{z_{m}}) = \\sigma^{2} = (1 - \\varepsilon)\\int^{\\infty}_{-\\infty}r^{2}\\phi(r)dr + (\\varepsilon)\\int^{\\infty}_{-\\infty}r^{2}\\phi(r/a)dr = 1 - \\varepsilon + \\varepsilon a^{2} $$\n",
    "\n",
    "Here $a^{2}=25$ so:\n",
    "\n",
    "$$  lim_{m\\rightarrow\\infty}V(\\sqrt{m}\\bar{z_{m}}) = 1 - \\varepsilon + 25\\varepsilon = 1 + 24\\varepsilon $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2\n",
    "Express the asymptotic variance of the normalized sample median at G as a function of $\\varepsilon$,\n",
    "which is given by:\n",
    "$$ lim_{m\\rightarrow\\infty}V(\\sqrt{m}\\theta_{m}) = \\frac{E[\\psi^{2}(r)]}{(E[\\psi'(r)])^{2}} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sample median has $\\psi_{med}(r)=sign(r)$ because $\\rho(r) = |r|$.\n",
    "\n",
    "With $E[\\psi^{2}_{med}(r)] = 1$ and $\\psi'_{med}(r) = 2\\partial(r)$:\n",
    "$$ E[$\\psi'_{med}(r)] = E[2\\partial(r)] = \\int^{\\infty}_{-\\infty}2\\partial(r)f(r)dr = 2f(0) $$\n",
    "\n",
    "$$ lim_{m\\rightarrow\\infty}V(\\sqrt{m}\\theta_{m}) = \\frac{1}{(2f(0))^{2}} = \\frac{1}{(4f(0))}$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ f(0) = \\frac{(1 - \\varepsilon + \\varepsilon/5)}{\\sqrt{2\\pi}} = (1-4\\varepsilon/5)(1/\\sqrt{2\\pi)} $$\n",
    "$$ 4f^{2}(0) = \\frac{2\\pi}{4-32\\varepsilon/5+64\\varepsilon^{2}/25} = \\frac{\\pi}{2-16\\varepsilon/5+32\\varepsilon^{2}/25}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3\n",
    "On the same figure, plot the asymptotic variance of the sample median and that of the sample mean versus $\\varepsilon$ for $\\varepsilon$ varying from 0 to 1. Find the value of $\\varepsilon$ at which the sample median becomes more efficient than the sample mean at G. Plot the density function of G at that value of $\\varepsilon$. Discuss the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgUAAAD8CAYAAADnntf1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VGW+x/HPk0YNmECAUEMLIQQpiYgoAmJBZVVEUUAE\nwQX7Ct6rWNZV790VXXGvrFgoinXBskhZFBQBBZHeuyACIYEQSkJLMpnn/nEGRAQSyCQnM/m+X695\nZeacMzPf8wrM/PKcpxhrLSIiIiIhbgcQERGR0kFFgYiIiAAqCkRERMRHRYGIiIgAKgpERETER0WB\niIiIACoKRERExEdFgYiIiAAqCkRERMQnrCTfrHr16jYuLq4k31JEJOAtW7Zsn7U2pgjPrxEWFjYO\nSEJ/DJZlXmCtx+O5Nzk5ee+ZDih0UWCMCQWWAqnW2u7GmGhgEhAHbAd6WWsPnOs14uLiWLp0aWHf\nUkREAGPML0V5flhY2LhatWo1j4mJORASEqK57csor9drMjIyEtPT08cBN53pmPOpGP8EbDjl8XBg\ntrW2KTDb91hEREqfpJiYmCwVBGVbSEiIjYmJOYTTYnTmYwrzQsaYusCNwLhTNt8MvOe7/x5wywXm\nFBGR4hWigkDAKQw4x3d/YVsK/g94HOd6xAk1rbVpvvvpQM0LSigiIiKlQoFFgTGmO7DXWrvsbMdY\nZ/3lM1ahxpjBxpilxpilGRkZF55URESkCIYPH16roGP27dsXOmLEiJOdOrdv3x7erVu3RsWbrPQo\nTEvB5cBNxpjtwETgKmPMh8AeY0wsgO/nGXsyWmvHWGtTrLUpMTEX3HlWRESkSEaNGhVb0DGZmZmh\n48ePr3HicVxcXN5XX321rXiTlR4FFgXW2iettXWttXHAncC31tq7gKlAf99h/YEpxZZSREQC2tVX\nX924RYsWzZs0adLilVdeqQ7g8Xjo2bNnXNOmTVvEx8cnPv/88zXWrVtXLjExsfmJ561Zs+bk4zp1\n6rR88MEH6yQkJCQmJSU1nz9/fsUrrriiab169ZJefvnlGIDp06dHpqSkNOvcuXOTuLi4pD59+tTP\nz8/ngQceqJOTkxOSkJCQeNNNNzUEeO6552o2bdq0RdOmTVu88MILNQAee+yxujt37iyXkJCQOGTI\nkLqbNm2KaNq0aYsTeQcPHlz3RN6//vWvNU4/z3bt2jUbNGhQvaSkpOaNGjVqMW/evIrXXntt4wYN\nGiQ98sgjtU8c98Ybb0S3bNmyeUJCQmKfPn0aeDweAPr27Vs/KSmpeZMmTVoMHTr05PF16tRpOXTo\n0NqJiYnN4+PjE1esWFG+GH5NRZqnYATwiTFmEPAL0Ms/kUREpFi1a9fsd9u6dTvICy/suaD9ixdv\nKugtP/roo+01a9bMP3z4sGnTpk3iXXfddWDLli3l0tLSwrds2bIOnKb76tWr50dGRub/8MMPFTp0\n6HDs7bffrt63b9/ME69Tv3793I0bN64fNGhQvYEDB8YtWrRo47Fjx0JatmzZ4vHHH88AWLNmTaUV\nK1asjY+Pz73yyiubvv/++1FvvPFG6oQJE2ps3LhxPcD3339f8eOPP662bNmyDdZakpOTm3ft2jV7\n5MiRu7p3717hxHGbNm2KOPHeI0eOjNmxY0fE+vXr14WHh7Nnz57QM51rRESEd+3atRv+53/+p8bt\nt9/eZMmSJRtq1KjhiYuLa/nUU0/t2b17d/hnn30WvXTp0o3lypWzd911V/233nqr2kMPPZT56quv\nptasWTPf4/HQoUOHZosWLapw6aWXHgOoXr26Z/369RtGjBgRM2LEiJqTJk0q0lDVMzmvSSystXOt\ntd199zOttV2ttU2ttVdba/f7O5yISKDLOJLBk988yebMzW5HcdVLL71Us1mzZonJycnN09PTw9et\nW1c+ISEhZ+fOneX69+9f77PPPqsSFRWVDzBgwIB9Y8eOre7xeJgyZUrUoEGDThYFvXr1OgjQsmXL\no23btj0SFRXlrV27ticiIsK7b9++UN++I4mJiblhYWH06tVr//fff1/59Dxz586tfMMNNxysUqWK\nt2rVqt4bb7zxwJw5cyLPdQ7ffvttlSFDhuwLDw8HoGbNmvlnOq5Hjx4HAVq1anWsSZMmxxo0aJBX\noUIFW69evZxt27ZFfPXVV5Fr166t2KpVq+YJCQmJ8+fPr7Jt27ZyAO+99150YmJi88TExMQtW7aU\nX7Vq1ckWgT59+hwAaNeu3dGdO3eWO69fQCGV6IyGIiJlzdLdSxmxYATXN72e+GrxbsdxFPSXfVH3\nn2b69OmR8+bNi1y6dOnGyMhIb7t27ZodO3YsJCYmJn/t2rXrJ0+eXOWtt96KmTRpUvSnn366vX//\n/gdeeuml2hMnTsxu2bLl0Vq1ap388i1fvrwFCAkJISIi4mQH95CQEPLy8gyAMeY373/64+J2asZy\n5cr9JqPH4zHWWnP77bdnjh49OvXU523cuDHi9ddfr7ls2bINMTEx+T179ow7fvx4yOmvGxYWZj0e\nT7GclKa7FBEpRsvTlgPQulZrl5O45+DBg6FVq1bNj4yM9K5YsaL8qlWrKgGkpaWF5efnM2DAgIMv\nvvhi6po1ayoCVKxY0Xbq1OnQsGHD6g8YMGDf+b7fmjVrKm3cuDEiPz+fzz77LLpjx47Z4HyZ5uTk\nGIAuXbocnjFjxkXZ2dkhWVlZITNmzIjq0qVLdtWqVfOPHDlyxu/Grl27Zr399tvV8/LyAM56+aAg\n3bp1y5o+fXpUampq2InX2bx5c8SBAwdCK1So4I2Ojs7fuXNn2Ny5c6teyOsXhVoKRESK0fL05TSN\nbkqVclXcjuKanj17HhozZkxMo0aNWjRq1Oh4q1atjoAz3G/QoEFxXq/XALzwwgu7Tjzn7rvv3v/V\nV19F3XrrrVnn+35JSUlH7rvvvvrbt28v36FDh6x+/fodBOjbt29G8+bNE5OSko5OnTr15z59+mS2\nbdu2OUC/fv0yLr/88mMAycnJh5s2bdriqquuOjRs2LCTI+uGDh2asXnz5nIJCQktwsLCbP/+/TOe\neuqp8x5rn5ycfPyZZ55J7dq1a7zX6yU8PNyOGjVqR9euXY8kJSUdbdy4cVJsbGxucnLy4fN97aIy\nzhQDJSMlJcVq7QMRKUsavtaQS+tcysTbJl7waxhjlllrUy70+atWrdreqlWr8/6L203PPvtszUOH\nDoW+9tpru8/nedOnT48cOXJkzTlz5vxUXNkC3apVq6q3atUq7kz71FIgIlJM9h/bz/aD27k/5X63\nowSUa665pvEvv/xSbt68eWW7d6YLVBSIiBSTFWkrAGgb29blJIHl66+/3nqhz+3evXt29+7ds/2Z\npyxRR0MRkWJyopNhm1ptXE4iUjgqCkREisny9OU0qNqAahWruR1FpFBUFIiIFJPlact16UACiooC\nEZFikJWTxZbMLbp0IAFFRYGISDFYtnsZFku7Ou3cjiI+gbZ08qhRo6rdfffd9QFefvnlmNdff73Y\nr0OpKBARKQaLUxcDcEmdS1xOIicE8tLJjz/+eMZDDz2UWfCRRaOiQESkGCxKXUST6CZEV4h2O0qp\noKWTC7d08muvvVYtLi4uqWXLls1/+OGHkws5DRs2rPazzz5bE2DkyJHVk5KSmjdr1izxuuuua5yd\nnR0C0LNnz7gBAwbUa9OmTULdunVbvvvuu1Hn+3vSPAUiIsVgcepiOsV1cjvG7wycMrDe2r1rK/rz\nNZNqJB195+Z3dp7rGC2dXPDSyX/4wx+yRowYUXvZsmUboqOj8zt06NAsKSnp6Omv37dv3wOPPfbY\nPoBHHnmk9qhRo6o//fTTewH27NkTvnTp0o0rV64s36NHjyb33HPPgfP5XaqlQETEz1KzUknNTqVd\nbfUnOEFLJxe8dPJ3331XqX379tm1a9f2lC9f3t566637z/T6y5Ytq5CcnNwsPj4+8fPPP6+2bt26\nk8sr33TTTQdDQ0NJTk4+npmZGV7Ar+V31FIgIuJnS3YvAeDSupe6nOT3CvqLvjho6eRfM55r6eQP\nPvjgosK8/uDBgxt+9tlnP1122WXHRo0aVW3evHkni5kT7w1wIWsbqaVARMTPFqcuJiwkrEwvl3wq\nLZ38W2dbOvnKK688smjRosj09PTQnJwcM3ny5DP2CTh69GhI/fr183JycszEiRP92mmlwKLAGFPe\nGLPYGLPKGLPOGPO8b/tzxphUY8xK3+0GfwYTEQlUi1MX06pmK8qHlS/44DKgZ8+ehzwej2nUqFGL\n//7v/65z6tLJV1xxRbOEhITEfv36NTp96WRjDEVZOrlx48ZJ9evXzzl96eSbbrqp4RVXXHH0xNLJ\nycnJzU8snVyrVq38E0snDxkypO6przt06NCMunXr5iYkJLRo1qxZ4vjx4y/oC/nUpZPj4+MTr7rq\nqvidO3eGN2jQIO+JJ57Y3b59++YpKSkJ8fHxx8/0/OHDh+9u165d85SUlISmTZue8ZgLVeDSycZp\nd6lkrT1sjAkH5gN/AroBh621rxT2zbR0sogEO6/1EvVSFH1b9uWNG9/wy2tq6eTC09LJBSvS0snW\nqRoO+x6G+27nf6FCRKQM2LRvE1k5WZq0qAi0dLJ7CtWnwBgTaoxZCewFvrbWLvLtetgYs9oY844x\n5rzHQ4qIBJsTkxapKLhwX3/99dbNmzevj42N9Zzvc7t3756tVoILV6iiwFqbb61tDdQF2hljkoA3\ngUZAayANGHmm5xpjBhtjlhpjlmZkZPgptohI6bRw10KqlKtCs2rN3I5yKq/X6y3ZLvhSKvn+HXjP\ntv+8Rh9Yaw8Cc4Bu1to9vmLBC4wFzlgWW2vHWGtTrLUpMTExZzpERCRozN8xnw71OhAackEd04vL\n2oyMjKoqDMo2r9drMjIyqgJrz3ZMgX0KjDExQJ619qAxpgJwDfCSMSbWWpvmO6zHud5ERKQsOHDs\nAOsy1nFn0p1uR/kNj8dzb3p6+rj09PQkNBS9LPMCaz0ez71nO6AwkxfFAu8ZY0Jx/jF9Yq2dboz5\nwBjTGqfT4XZgiB8Ci4gErIW7FgJwRf0rXE7yW8nJyXuBm9zOIaVfYUYfrAZ+tyC4tbZfsSQSEQlQ\n83fMJywkTJ0MJWCpGUlExE8W7FxAm1ptqBju1/WGREqMigIRET/Izc9lceriUnfpQOR8qCgQEfGD\n5WnLOe45zuX1Lnc7isgFU1EgIuIHC3YsAODy+ioKJHCpKBAR8YMFOxfQOKoxtSrXcjuKyAVTUSAi\nUkTWWubvmK9WAgl4KgpERIpofcZ6Mo5m0KlBJ7ejiBSJigIRkSKas30OAF3iuricRKRoVBSIiBTR\nnO1zaFC1AQ2jGrodRaRIVBSIiBSB13qZu30uXRqqlUACn4oCEZEiWLt3LfuP7delAwkKKgpERIpg\nzs9Of4LOcZ3dDSLiByoKRESKYM72OTSKakT9qvXdjiJSZCoKREQuUL43n3m/zNOlAwkaKgpERC7Q\nqj2rOHj8oIoCCRoqCkRELtDsbbMBNPJAgoaKAhGRCzRz60ySaiRRO7K221FE/KLAosAYU94Ys9gY\ns8oYs84Y87xve7Qx5mtjzBbfz6jijysiUjoczTvK9zu+57rG17kdRcRvCtNSkANcZa1tBbQGuhlj\n2gPDgdnW2qbAbN9jEZEyYd72eeTm53Jt42vdjiLiNwUWBdZx2Pcw3HezwM3Ae77t7wG3FEtCEZFS\naObWmZQPK0/H+h3djiLiN4XqU2CMCTXGrAT2Al9baxcBNa21ab5D0oGaxZRRRKTUmbV1Flc2uJIK\n4RXcjiLiN4UqCqy1+dba1kBdoJ0xJum0/Ran9eB3jDGDjTFLjTFLMzIyihxYRMRtOw/tZMO+DepP\nIEHnvEYfWGsPAnOAbsAeY0wsgO/n3rM8Z4y1NsVamxITE1PUvCIirpu1dRaA+hNI0CnM6IMYY8xF\nvvsVgGuAjcBUoL/vsP7AlOIKKSJSmszcOpM6kXVoEdPC7SgifhVWiGNigfeMMaE4RcQn1trpxpiF\nwCfGmEHAL0CvYswpIlIq5OXnMWvrLHo274kxxu04In5VYFFgrV0NtDnD9kyga3GEEhEprb7f8T2H\ncg7xh2Z/cDuKiN9pRkMRkfMwbdM0yoWW45pG17gdRcTvVBSIiBSStZZpm6dxVcOrqBRRye04In6n\nokBEpJA27tvI1gNbuanZTW5HESkWKgpERApp2uZpAHSP7+5yEpHioaJARKSQpm6aSptabahbpa7b\nUUSKhYoCEZFC2Hd0Hwt3LeQP8Rp1IMFLRYGISCHM2DIDr/VqKKIENRUFIiKFMHnjZGpH1qZtbFu3\no4gUGxUFIiIFyM7J5sstX3Jb89sIMfrYlOClf90iIgWYsWUGOfk53JZ4m9tRRIqVigIRkQJ8uv5T\nalWuRYd6HdyOIlKsVBSIiJzDkdwjzNgyg57NexIaEup2HJFipaJAROQcZmyZwTHPMV06kDJBRYGI\nyDl8tuEzalSqQcf6Hd2OIlLsVBSIiJzF0byj/Gfzf7g14VZdOpAyQUWBiMhZzNgygyN5R3TpQMoM\nFQUiImfx4eoPqVW5Fp3jOrsdRaREFFgUGGPqGWPmGGPWG2PWGWP+5Nv+nDEm1Riz0ne7ofjjioiU\njMyjmczYMoM+SX106UDKjLBCHOMBHrPWLjfGRALLjDFf+/b9w1r7SvHFExFxx6frPyXPm0e/Vv3c\njiJSYgosCqy1aUCa7362MWYDUKe4g4mIuOnD1R/SIqYFrWq2cjuKSIk5rz4Fxpg4oA2wyLfpYWPM\namPMO8aYKD9nExFxxbYD21iwcwF3XXwXxpiivdicObBggX+CiRSzQhcFxpjKwOfAo9baLOBNoBHQ\nGqclYeRZnjfYGLPUGLM0IyPDD5FFRIrXx2s+BqBPyz4X9gI7dsCBA8799HQYO9ZPyUSKV6GKAmNM\nOE5B8JG19t8A1to91tp8a60XGAu0O9NzrbVjrLUp1tqUmJgYf+UWESkW1lo+XP0hneM6U79q/cI/\n0euFL7+Em26Chg3hrbec7bffDuPHF09YET8rzOgDA4wHNlhrXz1le+wph/UA1vo/nohIyVq4ayGb\nMjfR7+JCdjD0euGll6BJE7jhBli0CIYPhz6+VoawMAjV6AUJDIUZfXA50A9YY4xZ6dv2FNDbGNMa\nsMB2YEixJBQRKUHjlo+jckRlerXodfaDrIUtWyA+HkJCYPp0qF8fXnwRevSAiIiSCyziR4UZfTAf\nOFNPmxn+jyMi4p6snCwmrZtE35Z9qRxR+QwHZMGHH8Kbb8LmzbBrF8TEwKxZUKFCyQcW8TPNaCgi\n4vOvNf/iaN5R/tj2j7/d8fPPcN99UKcOPPgghIfD669DpUrOfhUEEiQKc/lARKRMGLt8LBfXvJiU\n2ilw/LjTMlCjBmRnw3vvwR13wP33Q7t2UNShiiKlkFoKRESAFWkrWJa2jD/WvwXzxBNQty489piz\n8+KLYc8emDABLr1UBYEELbUUiIgA4774M+W8IfTt+QLkhcItt8DAgb8eUKWKe+FESoiKAhEpu/bu\nhZgYsnMP80H6LG7fWo6oJ4fDvfdC7dpupxMpcSoKRKRssdaZevjNN+GLL+C773gvdBnZIXk8/Ne5\n0KCD2wlFXKOiQETKhsOHYdy4X4cTRkfDo4/irVOb16feQ7s67WingkDKOHU0FJHgdmINgvx8ePpp\nqFYN3n8fUlPh73/nm7xNbMrcxMPtHnY3p0gpoJYCEQk+R47AxIlOq0B+PixfDlWrOrMQntZXYNSi\nUdSsVJPbE293KaxI6aGWAhEJHps3w5/+5EwydO+9zlwD997rrE8AvysItu7fyowtMxiSPIRyYeVc\nCCxSuqilQEQCW26u03mwXDn46iundeC225xJhq644pxzCvxz8T8JDQllSIqWbhEBtRSISKDasQOe\necZZiOiDD5xtAwc66xF8/DF07HjOgiDzaCZjl4+lT8s+1I7U8EMRUEuBiAQSa2HmTHjjDfjPf5zH\nN94ICQnO/sqVnVshvLn0TY7mHeW/LvuvYgwsElhUFIhI6ZeT41weAHj0UWdEwfDhMHgwNGhw3i93\nLO8YoxaN4oamN9CyZks/hxUJXCoKRKR0shYWLHD6CMyeDdu2QcWKMG2aUwhERFzwS09YOYGMoxk8\ncfkTfgwsEvhUFIhI6ZKVBR9+6BQDa9c6aw707w/HjjlFQdOmRXr5fG8+ryx8hUvrXErH+h39FFok\nOKgoEJHSweOBsDBYsgQefBDatoWxY6F3b6hUyW9v88m6T9h2YBsvX/0yRqsdivxGgaMPjDH1jDFz\njDHrjTHrjDF/8m2PNsZ8bYzZ4vsZVfxxRSSoHD/ujBzo0AH++7+dbVddBUuXOrd77/VrQZDvzeeF\n714gqUYSPZr38NvrigSLwgxJ9ACPWWsTgfbAg8aYRGA4MNta2xSY7XssIlKwn35yioC6deHuuyEz\nExITnX3GQHLyOYcTXqhJ6yaxcd9G/tLpL4QYjcgWOV2Blw+stWlAmu9+tjFmA1AHuBno7DvsPWAu\noF47InJmXi+E+L6IH38cpk6FHj2cSYa6dCmWIuBU+d58Xpj3Ai1rtOTW5rcW63uJBKrzKpWNMXFA\nG2ARUNNXMACkAzX9mkxEgsPu3fDCCxAX56w9APDyy87kQ59+6lwuKIFr+xPXTmRT5ia1EoicQ6E7\nGhpjKgOfA49aa7NO7aBjrbXGGHuW5w0GBgPUr1+/aGlFJDBYC99+64wg+OILZ1Gia66Bo0ed/U2a\nlGicE30JLq55sfoSiJxDoYoCY0w4TkHwkbX2377Ne4wxsdbaNGNMLLD3TM+11o4BxgCkpKScsXAQ\nkSBhrfNXf1oaXHedszLh0KEwZEiJFwKn+nD1h2zO3MznvT5XK4HIORRYFBinSWA8sMFa++opu6YC\n/YERvp9TiiWhiJR+S5Y4Uw8fOOC0DNSuDbNmwWWXQYUKrkY7lneMZ+Y8wyW1L6FHgloJRM6lMC0F\nlwP9gDXGmJW+bU/hFAOfGGMGAb8AvYonooiUSkeOwMSJziWCZcucoYP9+jmXCkJDnb4CpcCoRaPY\nlbWLD3t8qHkJRApQmNEH84Gz/U/q6t84IhIwXnkFnnsOkpJg9Gi46y5n9sFSZN/Rffxt/t/oHt+d\nTnGd3I4jUurp4pqIFCw3FyZNcoYOTp3qbBs8GL7/HlavhgceKHUFAcBfv/srh3MPM6LrCLejiAQE\nTXMsImf3yy8wZgyMHw979kDDhk6BABAb69xKqW0HtjF6yWgGth5Iixot3I4jEhBUFIjImeXnO9MP\np6fDjTc6kwxdd92vExCVcsNmDiMiNILnuzzvdhSRgKGiQEQcGRnwzjvw5ZfOUsWhofDee86qhA0a\nuJ3uvHy55UumbJrCiK4jqB1Z2+04IgEjMEp+ESke1jr9Avr0cdYhGO5bwmSvb9qRq68OuIIgx5PD\nI189Qny1eIZeNtTtOCIBRS0FImXZtGlw881OJ8EhQ+C++35dmChAvbrwVX7a/xMz75pJRGiE23FE\nAoqKApGyZNUqZ16B5s3hT3+Cbt2cSwa9evl1iWK37Dy0k//9/n/pkdCDaxtf63YckYCjywciwe74\ncfjgA6fTYOvWTj+BPXucfRERcM89QVEQWGt5cMaDWGt59bpXC36CiPyOWgpEgt3tt8P06RAfD6++\nCv37Q3S026n87tP1nzJt8zT+fs3fibsozu04IgFJLQUiwcTjcdYeuPHGXzsLPv44fPMNbNzoLE4U\nhAVB5tFMHv7yYZJjk3m0/aNuxxEJWGopEAkGu3fDuHEwdizs2gV16sCWLVCjBnTs6Ha6YvfYrMfY\nf2w/s+6aRViIPtZELpT+94gEup07oVEjp5Xg2mvhn/+E7t0hrGz89561dRbvrXqPp654ila1Wrkd\nRySglY1PDZFgcuAATJjgTDb0t79BvXpOX4Hrr4cmTdxOV6L2H9vPPVPuoXn15vy505/djiMS8FQU\niAQCa2HJEmc44cSJzoiCzp3B63WmHX74YbcTljhrLff/5372HtnLtN7TKB9W3u1IIgFPHQ1FAsHf\n/gaXXgqffuqMHli5EubMCZh1CIrDx2s+5pN1n/B85+dpG9vW7TgiQUEtBSKl0YYNTqtA375OMdCj\nB1x0EfTrVyqXKC5pOw7t4MEZD3J5vct54vIn3I4jEjRUFIiUFrm5znDCN96AefOciYWaN3eKgsTE\ngJ9+2F88Xg/9Jvcj3+bzfo/3CQ0JdTuSSNAosO3RGPOOMWavMWbtKdueM8akGmNW+m43FG9MkSCX\nn+986d9xB+zYASNGOKMK7r/f7WSlzrNznuW7X77jjRveoFFUI7fjiASVwrQUTABeB94/bfs/rLWv\n+D2RSFng9cLMmc7tH/9wlil+7DGIi4PrrivTfQXOZcaWGbw4/0XubXMv/Vr1czuOSNApsCiw1n5n\njIkr/igiZUBGhrMA0dtvw88/Q82azoyDtWurVaAAOw7toN/kfrSu1ZpR149yO45IUCrKnyMPG2NW\n+y4vRPktkUiwmjkT6taF4cOhQQOYNMm5VFC7ttvJSr0cTw69Pu1FXn4en97+KRXCK7gdSSQoXWhR\n8CbQCGgNpAEjz3agMWawMWapMWZpRkbGBb6dSADKynI6DX7xhfO4fXt44AFYv94ZTtirl9OZUM7J\nWssD/3mARamLePfmd2kSXbYmaBIpSRdUFFhr91hr8621XmAs0O4cx46x1qZYa1NiYmIuNKdI4Fi1\nCu67z2kBePBBmDzZ2V61qtN/oHlzd/MFmNcWvcY7K9/hmY7P0DOxp9txRILaBQ1JNMbEWmvTfA97\nAGvPdbxImdG/P7z/PpQvD717O/0ELrnE7VQBa9bWWTw26zF6JPTg+S7Pux1HJOgVWBQYY/4FdAaq\nG2N2AX8BOhtjWgMW2A4MKcaMIqXXTz/B+PHwzDNQqRJ06wZt2jjFQZS62hTFpn2buOOzO0iqkcT7\nPd4nxGhEhkhxK8zog95n2Dy+GLKIBAaPB6ZNc2Yc/PprZzhh165w9dVO64AUWfrhdLp91I3wkHCm\n3DmFyhFlNupgAAAW8klEQVSV3Y4kUiZoRkOR87FjB3ToAKmpUKcOPP883HuvRhD4UVZOFjd8dAN7\nj+xlTv85xF0U53YkkTJDRYHIuVgL334Lu3Y5lwTq1XMuEXTv7tzC9F/In3Lzc+n5SU9W71nNtN7T\naFfnrH2YRaQY6BNN5EwOHIAJE+Ctt2DzZmjSxFmMKCQExo1zO11Q8lovA6cM5Jtt3/Duze9yfdPr\n3Y4kUuao547I6d5807kcMGwYVKvmjCZYs0ZTDxcjay0P/udBPlrzEX+96q8MaD3A7UgiZZJaCkSO\nHIF//Qs6d3ZaBBITnUsF998PrVq5nS7oWWt59KtHeWvZWwy/fDhPXvGk25FEyiz96SNl1/r18Mgj\nTofBP/4RPv3U2d6pk3PZQAVBsbPW8vjXjzNq8SiGth/K37r+DWOM27FEyiy1FEjZk5/vrEQ4e7Yz\nzfBttznTD3fo4HayMsVay9PfPs0rC1/hwUseZOS1I1UQiLhMRYGUDb/8At98A4MGOfMKtGgB11wD\nAweCpt8ucV7rZehXQxm1eBSD2w5m1PWjVBCIlAIqCiR45ec7KxO++SbMmAHGwA03QGwsvPaa2+nK\nLI/Xwx+n/ZEJKycwtP1QtRCIlCLqUyDB6fvvnU6DN94IS5bAk0/Ctm1OQSCuyc3PpffnvZmwcgLP\ndXpOBYFIKaOWAgkO1sKCBVCunLMAUVwcNGoEL70Et9yiJYpLgaycLG7/9HZmbZ3FyGtHMuyyYW5H\nEpHTqCiQwJaVBR9+6FwiWLvWKQAmT3ZmHpw92+104rPz0E5u/PhG1mesZ/xN4xnYZqDbkUTkDFQU\nSOD685/hH/9w5hlo2xbGjtWCRKXQyvSV3PjxjWTnZDOj7wyubXyt25FE5CzUp0ACx/HjziRDHo/z\nuEIFuP12WLQIli51FiaqVMndjPIb/9n8Hzq+25FQE8qCgQtUEIiUcmopkNLvp5+cyYTefRf274eo\nKGdRoqeecjuZnIXXevnb93/j2TnP0rpWa6b3mU7tSK0kKVLaqSiQ0is11ZlHYNYsZ26BHj2cqYe7\ndHE7mZxDVk4Wd0++mymbptC3ZV/G/GEMFcMruh1LRApBRYGULrt3w9at0LGjM6nQ/v3w/PPOpYHa\n+kuztNuQsYEek3rw0/6f+L/r/o9HLn1EQw5FAkiBRYEx5h2gO7DXWpvk2xYNTALigO1AL2vtgeKL\nKUHNWvj2W2cEwRdfQN26zpwCERHOHAMSED5Y9QEPzHiAiuEVmX33bDrFdXI7koicp8J0NJwAdDtt\n23BgtrW2KTDb91jk/H3+OSQkwNVXw9y5znLFs2drmeIAkpWTRb/J/bj7i7tpG9uWZYOXqSAQCVAF\nthRYa78zxsSdtvlmoLPv/nvAXOAJP+aSYGWt89d/w4bO5YG8PKhe3RleeNttUL682wnlPCxJXULv\nz3vz88GfeaHzCzzV8SlCQ0LdjiUiF+hC/xyraa1N891PB2r6KY8EqyNHYNw4SEmBSy917gPccYcz\nE+Fdd6kgCCB5+Xm8MO8FOrzTgTxvHvMGzOPPnf6sgkAkwBW5o6G11hpj7Nn2G2MGA4MB6tevX9S3\nk0Dj8cBjj8F778GhQ5CUBKNHO0UAOIsUSUBZlb6KAVMGsDJ9Jb2TejP6htFEVYhyO5aI+MGFthTs\nMcbEAvh+7j3bgdbaMdbaFGttSoyWqC0bcnNh4ULnflgYrFvnLEz0/fewejU88ABUqeJuRjlvufm5\nPD/3eVLGppCWncbkOybzcc+PVRCIBJELbSmYCvQHRvh+TvFbIglcv/wCY8bA+PHOUMJdu6BGDWee\nAXUcDGgLdizg/v/cz5q9a+jbsi+vdXuNahWruR1LRPyswE9qY8y/gIVAM2PMLmPMIJxi4BpjzBbg\nat9jKatWr4Y//MFZlXDECGjXDqZMcToQggqCAJZxJIOBUwZyxbtXcPD4QabcOYUPb/1QBYFIkCrM\n6IOzrTDT1c9ZJJDs3QvHjkGDBk6/gCVL4MknYfBgUN+RgJfvzWfc8nE8OftJsnOzeeLyJ/jzlX+m\nUoTWlhAJZprRUArPWpg/35lk6LPPoFcvZ9nili2dSwVh+ucUDOZun8tjsx5jedpyOsd1ZvQNo0mM\nSXQ7loiUAH2KS+G8+y68+iqsXQtVqzprENx336/7VRAEvA0ZG3jimyeYtnka9arU46NbP6J3Um9N\nUyxShuiTXM5u3TpITHQuDyxaBOXKOfML3HmnligOInsO7+H5ec8zZtkYKkVUYkTXETxy6SNUCK/g\ndjQRKWEqCuS3jh+HTz5xLhH8+CP88ANcdhm89ppTFEjQyDiSwcsLXmb0ktHk5udyf8r9PNvpWWIq\naeiwSFmlokAcmZnw4ovOZYL9+yE+3rlckJDg7FdBEDT2Hd3HKz+8wuuLX+eY5xh9Wvbh2SufpWm1\npm5HExGXqSgoyzweZ6ni+vWdPgHjxsE11zj9Bbp00WyDQSb9cDqv/fgary95nSO5R7gz6U6e7fQs\nCdUT3I4mIqWEioKyaPdupwAYO9aZS2D5cqfz4K5dULmy2+nEzzbt28QrP7zC+6vfx+P1cFvibfyl\n0180okBEfkdFQVny44/wyivwxReQnw/XXuu0CpyggiCoLNy5kL//8He+2PgF5cLKMajNIIZdNowm\n0U3cjiYipZSKgmC3fz9UqODcliyBuXNh2DBnkqEm+nIINsfyjjFp3SRGLxnN0t1LiSofxTNXPsND\n7R6iRqUabscTkVJORUEwstYpAN58EyZOhH/+E+6917n98Y9aojgI/XzgZ95c+ibjV4xn/7H9NK/e\nnNevf53+rftTOUItQCJSOCoKgonHAxMmOMXA8uXO5YABA6BDB2d/BY07Dya5+blM3zydd1a8w4wt\nMwgxIdyScAsPXvIgneM6a9IhETlvKgqCwb59TofB0FAYORLCw+GNN6BvXy1RHIRWpa/i3ZXv8tGa\nj9h3dB+xlWN5uuPTDEkZQt0qdd2OJyIBTEVBoMrNhX//22kVWL3aGTlQqRLMmwcxMRpOGGT2HN7D\nJ+s+4d2V77IifQURoRHc3Oxm7ml9D9c0voawEP1XFpGi0ydJoElNdVoBxo1zVips2NBZndBaZ38N\ndSYLFplHM/n3hn8zad0k5myfg9d6aRvbln9e/096J/XW8sUi4ncqCgJBfr6zTHHlyrBpE4wYATfe\n6AwnvO46CAlxO6H4yaHjh5iyaQqT1k1i1tZZeLwemkY35emOT3NHiztoUaOF2xFFJIipKCjN9u6F\nd96Bt9+GW26Bf/zDmWlw+3aoV8/tdOInOw/tZOqmqUzdPJU5P88hz5tHg6oNGNZ+GHcm3UnrWq3V\naVBESoSKgtJowQIYPRo++wzy8qBzZ6cYAKevgAqCgGatZfWe1UzZNIUpm6awPG05APHV4nm0/aPc\n2vxWLq1zqQoBESlxRSoKjDHbgWwgH/BYa1P8EapMOnbs1yGDr74Ks2c7lwfuuw+aN3c3mxTZgWMH\n+GbbN8zaOouZW2eyM2snBsNl9S7jpatf4uZmN9OsejO3Y4pIGeePloIu1tp9fnidsmnlSmcEwccf\nw7JlzuqEo0bBRRc5owkkIHm8HhanLj5ZBCxOXYzXeqlaripdG3XlL53+Qvf47tSsXNPtqCIiJ+ny\ngRuOH4dPPnGKgR9/dGYY7N3bWakQoE4dd/PJecv35rMyfSXzfpnHvF/m8d0v33Hw+EFCTAjt6rTj\nmY7PcF2T62hXp52GD4pIqVXUTycLfGOMyQfettaO8UOm4JWbCxERkJkJAwc6aw/84x/Qvz9ERbmd\nTs5DXn4ey9OWnywC5u+YT1ZOFgBNoptwW/PbuLbxtXRt1JXoCtEupxURKZyiFgVXWGtTjTE1gK+N\nMRuttd+deoAxZjAwGKB+/fpFfLsA5PHAtGlOq4AxMHOm0xKwfDm0bKlJhgJExpEMFqUu4sddP568\nHck7AkBC9QR6J/WmU4NOdIrrRO3I2i6nFRG5MEUqCqy1qb6fe40xk4F2wHenHTMGGAOQkpJii/J+\nAWX3bhg71rmlpkLdujBkiDPJkDFw8cVuJ5SzyMvPY/We1SzctfBkAbD1wFYAQk0oF9e8mAGtB9Cp\nQSeubHCl+gWISNC44KLAGFMJCLHWZvvuXwu84LdkgcjrdX6GhMBbb8H//q8zudDo0c5kQ2G6llza\n5HhyWJexjuVpy1metpwV6StYmb6S457jANSqXIvL6l7G4OTBtK/bnuTYZCpFqAOoiASnonxL1QQm\n+8ZShwEfW2u/8kuqQLN/v7M64VtvwcsvOxMNPfSQ01egcWO304nPkdwjrNqz6jcFwNq9a/F4PQBE\nRkTSJrYN9yXfx2X1LqN93fbUq1JP8wWISJlxwUWBtXYb0MqPWQKLtbB4sdNXYNIkZ0TBZZdBZKSz\nv0YNrUPgkhxPDpsyN7Fu7zrW7l3L2oy1rNu7jm0HtmFxrmBVr1idtrFt+a/L/os2sW1oG9uWRlGN\nCDGaMlpEyi61Z5+vE30CPB6nRSA722kRuP9+aFV2ayQ3HPccZ+v+rWzYt8EpADLWsnbvWrZkbiHf\n5gNOH4Bm1ZvRNrYtd7e6m9a1WtM2ti11IuuoBUBE5DQqCgpr/Xrn8sDcubBiBYSHw5QpkJAAVaq4\nnS5oea2XnYd2silzE5szN//mtv3g9pN/+RsMjaMbk1QjiZ7Ne5JUI4mkGknEV4snIjTC5bMQEQkM\nKgrOJTcXJk92LhHMm+fMMXDbbXDoEERHQ7t2bicMCrn5uew4tINtB7bx84Gf+fngz2zZv4XNmZvZ\nkrmFnPyck8dGRkQSXy2e9nXbc3eru2lWrRnNqjcjoXoCFcMrungWIiKBT0XBmZy4RDBlCtx5JzRs\n6CxXPHAgxMS4nS7geK2X9MPpv/nS//ngzycfp2an4rXek8eHh4TTOLox8dXi6da4G/HV4mlWvRnx\n1eKpWammmv1FRIqJioIT8vOdiYXefBMuvxyGD4ebb4Yvv4Rrr3WGGcrvWGs5cPwAu7J2sfPQTnZm\n7fz1Z9bOk9tP/WsfoHZkbRpFNaJzXGcaXtSQRlGNaBjVkIYXNaR2ZG1CQ0JdOiMRkbJLRcHevTB+\nPIwZA9u3Q61acM01zr6ICOjWzdV4bvJ4Pew5vIf0w+mkHU4jLTvN+ZLP+u2X/9G8o795XqgJpXZk\nbepVrUfb2Lbc3Oxm50vf9+Xf4KIGlA8r79JZiYjI2ago6N0bvv0WunT5dY6B8HC3UxWrw7mHSctO\nO/lln344nbTstF/v+wqAfUf3nezId4LBEBsZS70q9UiqkcT1Ta6nXtV61KtSj3pV61G3Sl1qVa6l\nRX9ERAJQ2frkPnQIPvgA3nnHuSxQs6ZTCFSq5IwiCFBHco+QcTSDjCMZJ3/uPbLXuX/a9oyjGRzO\nPfy71wgLCaNW5VrEVo6lQdUGtK/TntjI2JPbTr0fHhrcRZOISFlVNoqCFSucvgIffwxHjkByMqSl\nOUVBcrLb6U7yWi9ZOVnsP7af/cf2c+DYgV/vHz9A5tHMM37RH/McO+PrRYRGEFMxhphKMcRUjKFJ\ndBNiKsb87ks+NjKW6ArRmrhHRKSMC/6iYMsWaNsWypd3LhXcfz9cckmxvZ3XesnOyeZQziEOHT/E\noZxDZ/2SP/3+weMHf9ML/3QVwyue/JKvUakGLWq0cB6f8sV/6s/IiEj11BcRkUILvqLgp5+cSYZy\ncuCf/4SmTZ0Wgm7dICrqnE/Nzc8lKyfr5Jf5oeOHnMenfMH/Zv/pj48fIjs3+5zvEWJCiCofRVSF\nKKIrRFOtYjWaRDchukI0UeWdbdEVok/uP7E9qkKUOueJiEixCtiiwOP1cDj3sHM7epDDs7/k8ORJ\nZK9ZxuEKIRzucAmHF77K4dwjHK5+mMM/zCc7N/vX55x2y87NPrky3rmUDytP1XJVqVq+KlXLVaVK\nuSrUql7L2eZ7fGLfiZ+nftFXKVdFzfQiIlIqBURR8OL3L/Lx2o9//QLPyf7duHcAWvpueIFFMGsR\n4HyRV46ofPIWGRFJ5YjK1Kpc6zfbz/SlXqVcld/c15S5IiISrAKiKIiqEEXT6Ka+L/TKVE7LpPKy\nNVS+8moqx7ekcmY2ken7qXx5FypXqPqbL/rKEZU1PE5ERKQQAuLb8r6U+7ivUS+YMMHpL7BlC1Sr\nBp0vh+Q73I4nIiISFAKiKCA3F5o1g337oEMHePZZZ2Gi8up4JyIi4i+BURRERMBrr0FSElx8sdtp\nREREglKRusEbY7oZYzYZY34yxgz3V6gz6tNHBYGIiEgxuuCiwBgTCowGrgcSgd7GmER/BRMREZGS\nVZSWgnbAT9babdbaXGAicLN/YomIiEhJK0pRUAfYecrjXb5tIiIiEoCKfWo9Y8xgY8xSY8zSjIyM\n4n47ERERuUBFKQpSgXqnPK7r2/Yb1tox1toUa21KTExMEd5OREREilNRioIlQFNjTENjTARwJzDV\nP7FERESkpF3wPAXWWo8x5iFgJhAKvGOtXee3ZCIiIlKiijR5kbV2BjDDT1lERETERcZaW3JvZkwG\n8MsFPr06sM+PcQKBzrls0DmXDUU55wbWWnXKkmJXokVBURhjllprU9zOUZJ0zmWDzrlsKIvnLIGn\n2IckioiISGBQUSAiIiJAYBUFY9wO4AKdc9mgcy4byuI5S4AJmD4FIiIiUrwCqaVAREREilGpKwqM\nMd2MMZuMMT8ZY4afYb8xxozy7V9tjGnrRk5/KsQ59/Wd6xpjzA/GmFZu5PSngs75lOMuMcZ4jDG3\nlWS+4lCYczbGdDbGrDTGrDPGzCvpjP5WiH/bVY0x04wxq3znfI8bOf3FGPOOMWavMWbtWfYH3eeX\nBBlrbam54cyMuBVoBEQAq4DE0465AfgSMEB7YJHbuUvgnDsAUb7715eFcz7luG9xJsi6ze3cJfB7\nvghYD9T3Pa7hdu4SOOengJd892OA/UCE29mLcM5XAm2BtWfZH1SfX7oF3620tRS0A36y1m6z1uYC\nE4GbTzvmZuB96/gRuMgYE1vSQf2owHO21v5grT3ge/gjzuJTgawwv2eAh4HPgb0lGa6YFOac+wD/\nttbuALDWBvp5F+acLRBpjDFAZZyiwFOyMf3HWvsdzjmcTbB9fkmQKW1FQR1g5ymPd/m2ne8xgeR8\nz2cQzl8agazAczbG1AF6AG+WYK7iVJjfczwQZYyZa4xZZoy5u8TSFY/CnPPrQHNgN7AG+JO11lsy\n8VwRbJ9fEmSKtPaBlCxjTBecouAKt7OUgP8DnrDWep0/IsuEMCAZ6ApUABYaY3601m52N1axug5Y\nCVwFNAa+NsZ8b63NcjeWSNlU2oqCVKDeKY/r+rad7zGBpFDnY4y5GBgHXG+tzSyhbMWlMOecAkz0\nFQTVgRuMMR5r7RclE9HvCnPOu4BMa+0R4Igx5jugFRCoRUFhzvkeYIS11gI/GWN+BhKAxSUTscQF\n2+eXBJnSdvlgCdDUGNPQGBMB3AlMPe2YqcDdvl687YFD1tq0kg7qRwWeszGmPvBvoF+Q/NVY4Dlb\naxtaa+OstXHAZ8ADAVwQQOH+bU8BrjDGhBljKgKXAhtKOKc/Feacd+C0jGCMqQk0A7aVaMqSFWyf\nXxJkSlVLgbXWY4x5CJiJ03P5HWvtOmPMfb79b+H0RL8B+Ak4ivOXRsAq5Dk/C1QD3vD95eyxAbyw\nSiHPOagU5pyttRuMMV8BqwEvMM5ae8ahbYGgkL/n/wEmGGPW4PTIf8JaG7CrJxpj/gV0BqobY3YB\nfwHCITg/vyT4aEZDERERAUrf5QMRERFxiYoCERERAVQUiIiIiI+KAhEREQFUFIiIiIiPigIREREB\nVBSIiIiIj4oCERERAeD/AUza8OnOAWB/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x81b3a58>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import math\n",
    "import random\n",
    "\n",
    "# Create a range and apply the derived asymptotic variance functions derived above.\n",
    "x = np.arange(0,1, 0.001)\n",
    "smean = 1 + 24*x\n",
    "smedian = 3.14159/(2 - ((16*x)/5) + ((32*(x**2)/25)))\n",
    "\n",
    "#plt.plot(t, t, 'r--', t, t**2, 'bs', t, t**3, 'g^')\n",
    "plot1 = plt.plot(x, smean, 'r--', label=\"asymptotic mean\")\n",
    "plt.plot(x, smedian, 'g-', label=\"asymptotic median\")\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Find the $\\varepsilon$ points where the two areas intersect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intersection found at epsilon = 0.026\n",
      "Intersection found at epsilon = 0.924\n"
     ]
    }
   ],
   "source": [
    "idx = np.argwhere(np.diff(np.sign(smean - smedian)) != 0).reshape(-1) + 0\n",
    "for i in idx:\n",
    "    print(\"Intersection found at epsilon = \" + str(i/1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the probability density of this function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1154e5c0>"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdYAAAD8CAYAAADKWoW+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9//HXJwlh34kgWQirEEC2gAqKCy6gFqhYi8Wq\ntUppq23tt63W7t+21v5qW7VarXVrq9b6tS6oKFbEFaWA7AIa9gSQsG8BspzfH2fGTEKWCUxmJjPv\n5+ORx7lz77l3PkHkc8+5555jzjlEREQkMlJiHYCIiEgiUWIVERGJICVWERGRCFJiFRERiSAlVhER\nkQhSYhUREYkgJVYREZEIUmIVERGJICVWERGRCEqLdQA16dKli8vNzY11GCIiTcaiRYt2OOcyTvAa\nJ6WlpT0EDEINr9pUACvKysquHzFixPaaKsRlYs3NzWXhwoWxDkNEpMkws40neo20tLSHunXrNiAj\nI2N3SkqK5rutQUVFhRUXF+dt27btIWBiTXV0RyIiIkGDMjIy9imp1i4lJcVlZGTsxbfqa64TxXhE\nRCS+pSip1i/wZ1Rr/lRiFRERiSAlVhERaVKeeOKJ9rfddlu3uurcc889nTds2NAsWjGFCiuxmtl4\nM1tjZgVmdmsd9UaaWZmZXd7Qc0VERMIxbdq0vbfffvu2uuo8/vjjXTZt2tSgxFpaWnpigQXUm1jN\nLBW4D5gA5AFXmlleLfV+C7zW0HNFREQA1qxZk96zZ8+BU6ZMyc3NzR00ceLEns8//3zb4cOH9+/R\no8eguXPntrrnnns6X3311TkA48aN633vvfd2Bvjd737XZeLEiT0fffTRjitWrGh19dVX9+rfv3/e\ngQMHLDMzc/DWrVvTAN5+++1Wo0aNOgXgu9/9bvfJkyf3HD58eP/LLrusZ1lZGV/72teyBg0aNKBf\nv355v/vd77o09HcI53WbUUCBc24dgJk9BUwCPqpW7ybg38DI4zhXRETiyZM2olGu+yW3qL4qmzdv\nbvGvf/1r3YgRIzaceuqpA5544onOCxcuXP3kk092+PWvf33yxIkT9wTrPvbYYxvHjBnTv0+fPkfu\nu+++bvPnz1/VtWvX8vvvv/+kO++8c/PYsWMP1fd9n3zySYv58+evbtOmjbvzzju7tG/fvnzFihWr\nSkpKbOTIkf0/97nP7evfv//RcH/FcLqCM4HNIZ8LA/s+Y2aZwOeB+xt6rkiT4hy8/DJ8+9swcSL8\n9a9w8GCsoxJJKJmZmUdGjRpVkpqaSr9+/UrOO++8fSkpKQwfPvxQYWFh89C62dnZZbfddtuWSy+9\n9JTbb799c9euXcsb+n3jx4/f06ZNGwfw+uuvt3v66ac79+/fP2/YsGEDdu/enfbRRx+1aMj1IjVB\nxF3ALc65CjM7rguY2XRgOkBOTk6EwhKJoLIy+MY3fDINevFF+NnP4I03oH//2MUmEmlhtCwbS3p6\n+mev/KSkpNCiRQsHkJqaSnl5+TFJZvny5S3bt29fVlRUVOsz1dTUVFdRUQFASUlJlUZl69atK4Lb\nzjn7/e9/v2nKlCn7jjf+cFqsRUB2yOeswL5Q+cBTZrYBuBz4s5lNDvNcAJxzDzrn8p1z+RkZJzQr\nl0jkOQdTp/qk2qIF/OQn8Je/wNChsHUrnH8+bNgQ6yhFks7cuXNbzZkzp/2iRYs+uvfee7utXr06\nHaBNmzble/fuTQ3Wy8rKOvree++1Anj66ac71na9Cy64YO/999+fceTIEQNYtmxZ83379jXoDZpw\nKi8A+ppZTzNLB6YCM0MrOOd6OudynXO5wDPAN5xzz4dzrkiT8PDD8O9/Q4cOMGcO/O//wvTp8N57\ncNZZUFQEkyf7Vq2IREVJSYnNmDEj96GHHtqQm5tbevvtt2++5pprcisqKrj66qt33HTTTT2Cg5d+\n+tOfbvnBD36QM2jQoAGpqam1ToJx88037+jfv//hwYMHD+jbt+/AG264oUdpaWmDumLNufon2TCz\ni/HdvanAI865X5vZDADn3APV6j4GvOSce6a2c+v7vvz8fKe5giVuFBbCwIGwbx88+SRceWXV43v3\nwpAhsHEj3HsvfPObsYlTkpqZLXLO5Z/INZYuXbphyJAhOyIVUyJbunRplyFDhuTWdCysZ6zOuVnA\nrGr7Hqil7rX1nSvSpPzgBz6pTpzou4Ora98e/vhHuOwy+PGP4YorQI8zRJKWZl4SqUtBAfzrX9Cs\nGfzpT1Db4LzJk+GCC2DPHvjtb6Mbo4jEFSVWkbrceSdUVMBVV0Fdo9XN4Pbb/fZDD8GBA9GJT0Ti\njhKrSG22bYPHHvNJ8/vfr79+fj6MGeOfuf7tb40enojEJyVWkdo88ggcOeKfrQ4YEN453/mOL+++\n27d0RSTpKLGK1MS5ylbn9Onhnzd5su8y/uQTmDu3cWITkbimxCpSk//+Fz7+GLp2hQsvDP+8tDS4\n+mq//c9/Nk5sItJoRo0adcrbb7/d6kSuocQqUpNga3XaNJ8sG+JLX/LlM8/4rmQRiYpILft2opRY\nRao7ehSeespvX3NNw88fMMBPGLF3L7z6amRjE0lwa9asSe/Vq9fAqVOn9ujTp8/AMWPG9D1w4IDN\nmzev5ZAhQ/r369cv74ILLuhdXFycCr6Fed1112UPGjRowK9+9auuU6ZMyZ02bVrOkCFD+mdlZQ1+\n6aWX2n7hC1/I7dWr18ApU6bkBr9n2rRpOYMGDRrQp0+fgTfffHP3SP4OkZqEXyRxvPUW7N7tZ1s6\n9dTju8aXvgRLl/qZmiZNimx8ItFgjbRsnKt/cv9Nmza1ePzxx9eNHj1648UXX9zr73//e8e77rqr\n2x//+MdNl1xyyYHvfOc73W+55ZbujzzyyGaAo0eP2ooVK1YBTJkyJXfv3r1pixcvXv3kk092mDp1\nap833nhj9YgRI0pOPfXUAfPmzWs5evTokj/84Q9FXbt2LS8rK2P06NGnzJ8/v+Vpp51WEolfUS1W\nkeqef96Xkycf/zWCMzS99BIcPnziMYkkkczMzCOjR48uARg2bNihtWvXNt+/f3/qJZdccgDghhtu\n2PnBBx+0Cda/8sord4Wef8kll+wJLjPXuXPn0tAl6NauXdsc4G9/+1unvLy8AXl5eXmffPJJi6VL\nlzZoabi6qMUqEso5mBlYJ+JEWpo5OX7lmyVLfAv4oosiE59ItITRsmwsocvGpaamuj179tS6HBxA\n27Ztq7zbFrrMXPUl6MrKymz16tXp9957b9dFixatysjIKJ8yZUru4cOHI9bQVItVJNSHH/pJ97t3\nhxEn2BN26aW+fOmlE49LJIm1b9++vF27duWvvvpqG4CHH3648xlnnHHc05vt3r07tWXLlhWdOnUq\n37x5c9qbb77ZPnLRKrGKVPXCC76cOBFSTvB/j9DEGsYqUiJSu0cffXT9LbfcktWvX7+8ZcuWtbzj\njju2HO+1zjjjjJJBgwYd6t2796Arrrii14gRIyI6B2lYy8ZFm5aNk5gZMcK3WmfNggkTTuxaFRXQ\nrRsUF8OKFX4wlEgj0bJx0VXXsnFqsYoE7dgBixdD8+Zwzjknfr2UFLj4Yr/94osnfj0RaRKUWEWC\n5szxXbZnngktW0bmmsHE+p//ROZ6IhL3wkqsZjbezNaYWYGZ3VrD8UlmtszMlpjZQjM7M+TYBjNb\nHjwWyeBFIur11315wQWRu+a55/ryvff02o00BRUVFRW1LDosQYE/o1pX2ag3sZpZKnAfMAHIA640\ns7xq1eYAQ5xzQ4HrgIeqHT/XOTf0RPv/RRqNc5WtyvPPj9x1MzL8LExHjsC8eZG7rkjjWFFcXNxe\nybV2FRUVVlxc3B5YUVudcN5jHQUUOOfWAZjZU8Ak4KNgBedc6Iiq1kD8jYgSqUtBAWzcCJ07w7Bh\nkb32uHF+FqY5c+C88yJ7bZEIKisru37btm0Pbdu2bRB6VFibCmBFWVnZ9bVVCCexZgKbQz4XAqdV\nr2Rmnwd+A5wEXBJyyAGvm1k58Bfn3IM1fYmZTQemA+Tk5IQRlkgEBbuBx4078ddsqjvvPPjDH+CN\nNyJ7XZEIGzFixHZgYqzjaOoi9i+Ic+4551x/YDLwy5BDZwa6iCcA3zSzsbWc/6BzLt85l5+RkRGp\nsETC8/bbvgw+E42ksWP9CjkLFsC+fZG/vojElXASaxGQHfI5K7CvRs65t4FeZtYl8LkoUG4HnsN3\nLYvED+cqE+vYGu/7TkzbtjBqFJSXwzvvRP76IhJXwkmsC4C+ZtbTzNKBqcDM0Apm1sfMLLA9HGgO\n7DSz1mbWNrC/NXAhdTzwFYmJdetgyxbo0sUv+dYYgglbiVUk4dX7jNU5V2ZmNwKzgVTgEefcSjOb\nETj+ADAFuNrMSoES4IvOOWdmXYHnAjk3DXjSOacFKiW+BFurZ50F1kiDIc8MvIH27ruNc30RiRth\nrW7jnJsFzKq274GQ7d8Cv63hvHXAkBOMUaRxBVuRZ53VeN8xerRP2gsW+PdZW0RshSoRiTMaTi3S\nmM9Xgzp2hEGD4OhRn1xFJGEpsUpy27IF1q71A4yGNHLnirqDRZKCEqskt/ff9+Xpp/tXYhpTsKtZ\nA5hEEpoSqyS3YGI944zG/65gi3XePL+knIgkJCVWSW7RTKzZ2ZCVBXv3wurVjf99IhITSqySvI4e\nhUWL/PZpx8zS2ThOP92X8+dH5/tEJOqUWCV5LV7sV50ZMMCP2o2GYGL94IPofJ+IRJ0SqySvaHYD\nBwVbxkqsIglLiVWSVywS6/DhfvTxihWwf3/0vldEokaJVZJX8DlnsHs2Glq18u/LVlTAwoXR+14R\niRolVklO27f7hc1bt268ifdrowFMIglNiVWSU3Bawfx8SE1t+PmHt8OhWldPrJsGMIkktEaeakYk\nTgUT68iR4Z9zcBNsfAo2Pwc7A0mxTS846RzIngKZF4d3neB3qitYJCEpsUpy+u9/fTlqVHj1C1+E\n96ZC+SH/OaU5pDaHA+v8z7pHoN+3YPidkNKs7mv17Qvt2kFREWzdCieffPy/h4jEHXUFS/JxrmEt\n1o/vg3cm+6Ta/RI48xmYsgOm7ILxC2Hwz30y/fgemHsRHN5R9/VSUmDECL+tlW5EEo4SqySfjRth\nxw7o0gV69Ki9nnOw5FZYeCO4Cp9Az34RcqZAszaQkgqdRsDgn8G4udCiK3w6F147vf7kqu5gkYQV\nVmI1s/FmtsbMCszs1hqOTzKzZWa2xMwWmtmZ4Z4rEnWh3cBmtdfb8Dh89FuwNDj9UZ9Aa6ufMca3\nXjsOhQNr4YNrfTKuTTCxqsUqknDqTaxmlgrcB0wA8oArzSyvWrU5wBDn3FDgOuChBpwrEl3BVmJd\n3cAH1sGCb/rtUQ9Ar2vrv26rLBj7AqR3hC0vw+o/1F43P9+XCxb4lrGIJIxwWqyjgALn3Drn3FHg\nKWBSaAXn3AHnPvvXoTXgwj1XJOqCiTX4nLO6ijKYdxWU7Yfsy6HXdeFfu3UOnP6Y317yQ9hRyys1\nPXr4ruidO2HDhvCvLyJxL5zEmglsDvlcGNhXhZl93sxWAy/jW61hnxs4f3qgG3lhcXFxOLGLNJxz\n8OGHfru2xLril7Djfd8CHfWXuruLa5I1EU65GVwZvPtFKK1h6kIzPWcVSVARG7zknHvOOdcfmAz8\n8jjOf9A5l++cy8/IyIhUWCJVrV3r10Pt1g26dz/2+J6VsPJXgMEZ/4DmnY7ve4be4Qc2HdoEq35f\nc51gYg8uXSciCSGcxFoEZId8zgrsq5Fz7m2gl5l1aei5Io0umMRqa60u+7EfdNTna9D1nOP/ntR0\nGH6X3159J5RsPbZOMIZgC1pEEkI4iXUB0NfMeppZOjAVmBlawcz6mPn+MjMbDjQHdoZzrkhU1ZVY\nd3wAhc9DaisY/NMT/66TzoSsyVB2EJb/4tjjw4dXxqQBTCIJo97E6pwrA24EZgOrgKedcyvNbIaZ\nzQhUmwKsMLMl+FHAX3Rejec2xi8iEpZgYg2Oyg0KvrMK0P870DJCsyEN+Q1YKqx9CPauqnosO9sP\nYNq1CzZtisz3iUjMmYvDO+X8/Hy3UAM6JNKcg06dYM8eP51g6DPWLbPhzfH+VZmJ6yC9Q+S+979f\nh4IHIGsSjH2+6rGLLoLXXoN//xsuuyxy3ylJx8wWOefy668pjU0zL0nyWLfOJ9XqA5dcBSz9od/O\n+2Fkkyr4iSXSWkPhC7Cz2oQQes4qknCUWCV5BJNX8Nlm0JZXYfdiaNkd+t0Y+e9t2Q36fsNvr/5j\n1WOhz1lFJCEosUryWLzYl9UT65q7fXnKtyGtZeN8d7+b/LPWTf8HB0Ne7Q595SYOH8uISMMpsUry\nCLZYhw2r3Ld3FWx7DVJbQu/rG++7W2dDzhf8pBEf31u5PzcXOnaE4mL/3FdEmjwlVkkOoTMuhSbW\nNff4sufVxz8ZRLhOudmXBQ9C6QG/bQZDh/rtYItaRJo0JVZJDlu2+FZhhw6+lQhwdDes/7vfPuVb\njR9Dl1HQZTSU7oH1f6vcH0z0SqwiCUGJVZJDMGkNG1Y592/BQ37x8m4XQPsoLbrUP9BqXX1X5bJy\nSqwiCUWJVZJD9RHBFSHPOk/5TvTiyJoMrXPhQAFsne33KbGKJBQlVkkOoS1WgK2v+gny2/aF7uOj\nF0dKGvS5wW+ve9SXp5wCLVrAxo1+FiYRadKUWCU5VB+4tO4xX/a+HizK/xv0vNp/Z+ELcGQnpKXB\nkCH+2JIl0Y1FRCJOiVUSX3Au3pYtfevw8A4omumTW+5V0Y+nVRZ0uwgqjsKGJ/0+dQeLJAwlVkl8\nwVbgqadCaips/CdUlPrk1qqGNVmjofdXfLnuEV8qsYokDCVWSXzBxBp8XzTYDdzr2lhE42VOhPRO\nsHsJ7FqsxCqSQJRYJfGFDlzavQx2f+hXscmaGLuYUptD7jS/ve5RGDTIt6ZXr4aSktjFJSInTIlV\nEl9oizXYWu1xJaS2iFlIAPS+zpcbHod0g/79oaICli+PbVwickLCSqxmNt7M1phZgZndWsPxaWa2\nzMyWm9k8MxsScmxDYP8SM9MiqxJdJSWwahWkpMDA/j6JQWy7gYM6DoWOw/wMUEUvV3ZVa2SwSJNW\nb2I1s1TgPmACkAdcaWbVp6lZD5ztnBsM/BJ4sNrxc51zQ7UIr0TdypVQXu5HA+99B44U+1mWOsXJ\nX8XgqOSNT+k5q0iCCKfFOgoocM6tc84dBZ4CJoVWcM7Nc87tDnz8AMiKbJgixyn0+erGf/nt3GmV\n0xrGWo8rfLnlJRjUz2+rxSrSpIWTWDOBkAUkKQzsq81XgVdCPjvgdTNbZGbTGx6iyAn47FWbgVD4\nvN/O+WLs4qmuVRZknAXlh+GkwLJxy5b5VraINEkRHbxkZufiE+stIbvPdM4NxXclf9PMxtZy7nQz\nW2hmC4uLiyMZliSzYGLNPgJlB3wXcNvesY2puh5Tfbn3JcjJgUOH4JNPYhuTiBy3cBJrEZAd8jkr\nsK8KMzsVeAiY5JzbGdzvnCsKlNuB5/Bdy8dwzj3onMt3zuVnZGSE/xuI1Ka8HJYu9dvtAgm2Rxy1\nVoNyLgdL9ZPynxoYvqDnrCJNVjiJdQHQ18x6mlk6MBWYGVrBzHKAZ4EvO+c+Dtnf2szaBreBC4EV\nkQpepE5r18LBg5CVCQf/4/flXBHbmGrS4iToOg5cGfRs5vfpOatIk5VWXwXnXJmZ3QjMBlKBR5xz\nK81sRuD4A8BPgc7An80PCikLjADuCjwX2JcGPOmce7VRfhOR6oLJqV9XKC/yi4y3zoltTLXpMRW2\nvQZdNvjPSqwiTVa9iRXAOTcLmFVt3wMh29cD19dw3jpgSPX9IlER7E7NOuDLeOwGDsr+PCz4GrQP\nTA6xeDE4Fz+jl0UkbJp5SRJXsNXXeR1gkH15TMOpU3oHOHkCdAHatYTiYti6NdZRichxUGKVxBVs\nseaUwUljY7eSTbhyLgcDejb3n9UdLNIkKbFKYtq2DT79FFqn+VZgPLdWgzIvBUuD7nv9ZyVWkSZJ\niVUS02fvr1b4v+XZk2MaTljSO0K3cZDj/GclVpEmSYlVElOwG7hHBXQ+zc9w1BRkXwY9Att6l1Wk\nSVJilcQUbO31wCerpiJzkp8wNA0oKIB9+2IdkYg0kBKrJKbFH/qyB5D1+ZiG0iAtu0L3sZXLWCxb\nFtNwRKThlFgl8ezfDwVr/XQmAwdCu76xjqhhstQdLNKUKbFK4lm2zE+ukAX0bAKjgasLfc66aEFM\nQxGRhlNilcQT7AbOpWk9Xw1qnQ2D+/vtRe/GNhYRaTAlVkk87wemo+7bGToMjm0sx2vsVD9ZxOqN\ncPRorKMRkQZQYpXEs3ihL0+/oOnOtXvKFX4Ji7IKWKEBTCJNiRKrJJajR6Gg2G+ffU1sYzkR7fpD\n7zZ+++1nYhuLiDSIEqsklvkzodRBtxTodUGsozl+ZjBsmN+ePzu2sYhIgyixSmJ560lf5uVASmps\nYzlRZwamYVy2yo9yFpEmQYlVEsvC93w58qzYxhEJY6f6ct0R2K3nrCJNRViJ1czGm9kaMysws1tr\nOD7NzJaZ2XIzm2dmQ8I9VyRiDhXC6u1+e+yU2MYSCSd3h4xWcBiY92isoxGRMNWbWM0sFbgPmADk\nAVeaWV61auuBs51zg4FfAg824FyRyNj0PGwMbOefEdNQImbIQF++NzO2cYhI2MJpsY4CCpxz65xz\nR4GngEmhFZxz85xzuwMfP6ByptN6zxWJmPlP+dZd1w5w0kmxjiYyThvny5Xr4eDm2MYiImEJJ7Fm\nAqH/RxcG9tXmq8Arx3muyPE5uhc+/MBvjxgV21giKf80X64HCl+IaSgiEp6IDl4ys3PxifWW4zh3\nupktNLOFxcXFkQxLksHWV2Fdud8ecVpsY4mk4cN9uQElVpEmIpzEWgRkh3zOCuyrwsxOBR4CJjnn\ndjbkXADn3IPOuXznXH5GRkY4sYtUKpzpkw9Uvv+ZCLKzoVNHOACsftO3zEUkroWTWBcAfc2sp5ml\nA1OBKiMpzCwHeBb4snPu44acK3LCKkqh6OXKxDpiRCyjiSwzGB74fdaW+Za5iMS1ehOrc64MuBGY\nDawCnnbOrTSzGWY2I1Dtp0Bn4M9mtsTMFtZ1biP8HpLMtr8NW/f6Vl2XLr6Vl0iqdAfrvlQk3qWF\nU8k5NwuYVW3fAyHb1wPXh3uuSESFdgOPGNF0J96vTTCxrge2zPIt9JRmMQ1JRGqnmZekaXMOil7w\nSQcqk1AiCXZtb0yF0j2w/Z3YxiMidVJilaZtz3I4uBE2pvvPifR8Nah3b2jfHnaXw26gSN3BIvFM\niVWatsIXwAEbA3+VEzGxmlXtDi58QZPyi8QxJVZp2opmwi5g92Ho2BF69Ih1RI0jeMOwuRUc3AB7\nV8Q0HBGpnRKrNF2HimDXwqrdwIk2cCkomFiLOvtSk0WIxC0lVmm6gs8aPw20UhOxGzgo+Lt9ctCX\nSqwicUuJVZquYHLZ2NyXI0fGLpbG1rs3tGsH23fBvpa+pX6oMNZRiUgNlFilaSrdB5++ARh8tMXv\ny8+PaUiNKiWlstW6+1RfarIIkbikxCpN05ZX/EQJZfmwcxdkZEBOTqyjalzBG4ctXX2p7mCRuKTE\nKk1TMKns7O/L/PzEHbgUFOzq/vgAWApsn6tJ+UXikBKrND3lR2HLy347GZ6vBgV/x4VLoMsY32Lf\n8krd54hI1CmxStOz/S3/jLX9QFi21u9L5OerQT16+C7vXbvAneX3Fak7WCTeKLFK0xPsBu4+CRYt\n8tvJkFjNKlutweesW2b5FryIxA0lVmlagpPuAxweAvv2QWYmnHxybOOKlmBiXbHJt9hL98H2N2Ma\nkohUpcQqTcvuD/37my0zYfUBv++002IbUzQFE+uCBZA12W8XPh+7eETkGEqs0rRsDiSRrImwYKHf\nTsbEumgRdJ/otwufB1cRu5hEpIqwEquZjTezNWZWYGa31nC8v5m9b2ZHzOx71Y5tMLPlZrbEzBZG\nKnBJUoXP+jL7Mpg/328nU2I96STo2RMOHoQtzaFVDpRshR3zYx2ZiATUm1jNLBW4D5gA5AFXmlle\ntWq7gG8Bd9ZymXOdc0Odc0kwwkQazd7VsPcjSO8IbUfBsmVVZyRKFsEbifnz/Q0GVN5wiEjMhdNi\nHQUUOOfWOeeOAk8Bk0IrOOe2O+cWAKWNEKOIF0wemRNhyXIoK4NBg6BNm9jGFW2nn+7L0MS6+Vmt\n0SoSJ8JJrJnA5pDPhYF94XLA62a2yMymNyQ4kSo219ANPGpU7OKJlWCL9YMPoMtoaHESHFgHe5bH\nNi4RAaIzeOlM59xQfFfyN81sbE2VzGy6mS00s4XFxcVRCEualIMbYdciSGsN3S5IzuerQcOGQXo6\nrFoF+w9AZqADabO6g0XiQTiJtQjIDvmcFdgXFudcUaDcDjyH71quqd6Dzrl851x+RkZGuJeXZLH5\nOV92vxjSWvrWGiRnYm3eHIYO9V2/CxboOatInAknsS4A+ppZTzNLB6YCYa1XZWatzaxtcBu4EFhx\nvMFKEgu2xrIug8JC2LTJr086cGBs44qV4HPWDz6ArudBs3a+K3h/QWzjEpH6E6tzrgy4EZgNrAKe\nds6tNLMZZjYDwMy6mVkh8F3gx2ZWaGbtgK7Au2a2FPgv8LJz7tXG+mUkQZV8CsXvQko6ZF4C77/v\n959+uh8VnIyCifX99yE1HTI/5z+rO1gk5tLCqeScmwXMqrbvgZDtbfgu4ur2AUNOJEARP7OQg24X\nQrO2MG+e3z96dEzDiqkxY3w5bx5UVED2FNjwBGx6BvJ+ENvYRJJckt7uS5Oy6f98GXyWGGyxJnNi\nzcmBrCzYs8cPYjp5PKS1gV0L4MD6WEcnktSUWCW+lXzqF/ROaQbZk6GkBD780K/0kowDl0IFbyzm\nzfMDuoLdwcEbERGJCSVWiW+Fz/p5cLtd6GdcWrQISkth8GA/eCmZBbuD33vPlzlX+HLT07GJR0QA\nJVaJdxsDSaLHF30ZTCJnnBGbeOJJ9cTaPdgdvAj2r41dXCJJTolV4lfJVtj+VmA0cGAll3fe8eVZ\nZ8UurnjnK1y/AAAYbklEQVQxZAi0bg0FBfDpp5DaArICk0WoO1gkZpRYJX5t+jfgoPsESG8P5eXw\n7rv+mBIrpKVVPmf+rDv4C75Ud7BIzCixSvza9C9fBp8drlgBe/f6EbE5ObGLK54EbzDeftuXJ18E\naW1h92JNFiESI0qsEp8OFflJIVJbVI52DXYDj61xuunkdPbZvgwm1irdwWq1isSCEqvEp+Azwu4X\n+0khoDJ5qBu40mmnQbNmsGSJf6cVKgd6bXhCS8mJxIASq8SnDY/7ssdUXzqnFmtNWrXyS+c5V/mc\n9eSLoHlnvyj8nqWxjU8kCSmxSvzZu8q/MtKsHXS/1O8rKIBt2yAjA045JbbxxZvgjcZbb/kypRnk\nBFqt6x+PTUwiSUyJVeJPsLWa8wU/oxDA3Lm+POccP+uSVKr+nBUg9ypfbnwSKsqjH5NIElNilfji\nKvyzQYDcL1fuf+MNX557bvRjinejR0NqKixcCPv3+31dToc2vQLvAs+NbXwiSUaJVeJL8btwcCO0\nyoGTAoOUnKtMrOedF7vY4lXbtjBypH/PN9hqNatstao7WCSqlFglvqz/hy9zp4EF/nquXAnFxdC9\nO/TrF7vY4tn55/vy9dcr9+VO8+XmZ6HsUPRjEklSSqwSP8oPV75m0/Oqyv2hrVU9X63ZuHG+nDOn\ncl+7ftBpJJTth6IXYxOXSBIKK7Ga2XgzW2NmBWZ2aw3H+5vZ+2Z2xMy+15BzRT5T9BKU7oWOw6F9\nXuX+4MAldQPX7owzoGVLWL4ctm+v3B+8QVn3WEzCEklG9SZWM0sF7gMmAHnAlWaWV63aLuBbwJ3H\nca6IF/zHv2fIoKXSUj1fDUfz5nDmmX47+OcF0ONLfhGDrbPh4ObYxCaSZMJpsY4CCpxz65xzR4Gn\ngEmhFZxz251zC4DShp4rAsChQtj6in8HM/hsEGD+fNi3z7+72qNH7OJrCoLPWf/zn8p9LbpA1mTA\nwbpHYxKWSLIJJ7FmAqG3uoWBfeEI+1wzm25mC81sYXFxcZiXl4Sx9hH/qk3W56FFRuX+2bN9edFF\nsYmrKbngAl/Onl11KsPe1/tyXeDPWEQaVdwMXnLOPeicy3fO5WdkZNR/giSOinJY+7Df7nND1WNK\nrOEbMgS6doWiIj+SOqjbOGjdw7/GtG1O7eeLSESEk1iLgOyQz1mBfeE4kXMlWWz7DxzaBK17QteQ\n56g7dvhJD9LTK2cXktqlpMD48X77lVcq91sK9Pqq3177UPTjEkky4STWBUBfM+tpZunAVGBmmNc/\nkXMlWaz9qy/7XF/57ir4Z4XO+dVsWreOTWxNzYQJvgxNrAC9rvV/toXPweEdUQ9LJJnUm1idc2XA\njcBsYBXwtHNupZnNMLMZAGbWzcwKge8CPzazQjNrV9u5jfXLSBNU8ikUzgRLhZ7XVj02a5Yvg60w\nqd8FF/iW67vvVk5vCNA6G04eDxWlsOEfsYtPJAmE9YzVOTfLOdfPOdfbOffrwL4HnHMPBLa3Oeey\nnHPtnHMdAtv7ajtX5DPr/wauDLpfAq26V+4vK6tMrJdeGpvYmqJOnfwaraWlVSeLgMpBTJ88oEFM\nIo0obgYvSRKqKINP/uy3+0yveuz992HXLujbV8vENdQll/hyZrWnLpmfg1bZsP9j/16riDQKJVaJ\nncIX/EjVNn2g+4Sqx4JJ4XOf0zSGDTUp8Kr4Sy/5ifmDUtKg3zf99pp7oh+XSJJQYpXYWXO3L0/5\ndtVBSwAvBua2/dznohtTIhg4EHr18gsXvP9+1WO9b4DUlrD1Vdi3JjbxiSQ4JVaJjV2LoPgdaNYO\nel1T9dgnn8CaNdChA4wZE5v4mjIzmDzZb7/wQtVjzTtVLif38b3RjUskSSixSmysDrRWe18PzdpW\nPfbMM7689FJo1iy6cSWKYHfw889XnYUJ4JRv+XLdY3B0b1TDEkkGSqwSfSVbYdNTvvu3343HHv+/\nwNJxX/hCdONKJKNHQ5cuUFAAK1ZUPdZhkJ+Io+yAn+ZQRCJKiVWi75P7/fuUWZOhTc+qx9auhcWL\noW1buPDC2MSXCNLS4LLL/Pa//nXs8VO+7cvVd/n/FiISMUqsEl2l+2DNn/z2KTcfezzYDTxxIrRo\nEb24EtEXv+jLp546tjs481K/5u2hTbD+8ejHJpLAlFgluj6+F0r3wElj4aQzjz0e7Aa+/PLoxpWI\nzj7bT8q/di18+GHVY5YCebf57ZW3+4UQRCQilFglekoPwOo/+O1BPz32+KpVsGgRtGun1WwiITW1\n8galpu7gHl+ENr3hQAFsejq6sYkkMCVWiZ5P7ocjO6HLGVVXsQn6R2AO2yuugJYtoxtbopo61Zf/\n/GfVySLATxiRd6vfXvlrTXMoEiFKrBIdZYdg9Z1+e9BPj51NqaICHg886/vyl6MbWyIbM8ZPFlFY\nCG+8cezxnldDqyzYu9LPhCUiJ0yJVaKj4EE4vB065cPJNXTzvvUWbN4MublwZg3PXuX4mMG11/rt\nRx899nhqOgz4gd9e/nM9axWJACVWaXxH9/quRoBBP6l57t/gP/pXXeWXPZPIufpq/2f+3HOwZ8+x\nx3tf7yfn37MMNmiEsMiJ0r9g0vg++g0c2QEZZ/oVVqrbtQueDgye+cpXohtbMujRA847Dw4f9q/e\nVJfWEk79ld9e9mMoK4lufCIJJqzEambjzWyNmRWY2a01HDczuydwfJmZDQ85tsHMlpvZEjNbGMng\npQk4sMFPQgAw7Pc1t1YfewyOHPEjgXv1imZ0yeOrX/Xln/987DutAD2vgo5D4VAhrLkrurGJJJh6\nE6uZpQL3AROAPOBKM8urVm0C0DfwMx24v9rxc51zQ51z+ScesjQpS38IFUegx5egy6hjjzsHDzzg\nt2fMiG5syWTKFP9O6/Ll8M47xx63FBj2O7+98jdwuDi68YkkkHBarKOAAufcOufcUeApYFK1OpOA\nvzvvA6CDmZ0c4VilqdkxHzY+BSnNYejtNdeZM8evZpOZ6Sfdl8aRng7TA4vJ33dfzXW6nQ8nj4ey\n/bD8F9GLTSTBhJNYM4HNIZ8LA/vCreOA181skZlNP95ApYmpKIdFgVVU+t8MrXvUXO/OwCs4M2b4\n+W2l8Xzta37SiGef9a/f1GTY//Ot14L7YeeC6MYnkiCiMXjpTOfcUHx38TfNbGxNlcxsupktNLOF\nxcXqhmryPv4T7PwvtMyEgT+suc6yZTB7NrRqBV//enTjS0aZmX4mprIy+OMfa67TYbCfw9lVwPzr\nNUG/yHEIJ7EWAdkhn7MC+8Kq45wLltuB5/Bdy8dwzj3onMt3zuVnZGSEF73EpwPrYemP/PbI+/1i\n5jX5Q2B6w69+FTp3jk5sye6WW3z5l7/Azp011zn1F9A6179+s+r3UQtNJFGEk1gXAH3NrKeZpQNT\ngZnV6swErg6MDj4d2Ouc22pmrc2sLYCZtQYuBKotDikJxTn473QoPwQ9pkJWDa/XAKxfD0884d9Z\nvbmGVW6kcQwb5kdfHzwI995bc5201jDqL357xS9gf0H04hNJAPUmVudcGXAjMBtYBTztnFtpZjPM\nLDiMcxawDigA/gp8I7C/K/CumS0F/gu87Jx7NcK/g8ST9X+Dba9D884w4u7a6/3v//ouyS9/GXr2\nrL2eRN6tgTfm7roLdu+uuc7JF0LuVVB+ONAlrBmZRMJlrqZ32mIsPz/fLVyoV16bnH0fw6v5flTp\nGf/w70bWZM0ayMvzrdU1a/TuarQ5B+PGwdy5Psn+5jc11zu8A2YN9FNRDv45DP5ZVMOUhjGzRXql\nMT5o5iWJjLISePcLPqnmXAG502qv+5Of+En3r7tOSTUWzOCOO/z23XfDli0112vRBUY/AZh//ebT\nuVELUaQpU2KVyFh0kx/s0rYvnPbXmmdYAnjzTb+YecuW8OMfRzVECTFqFFx2GZSUwI9+VHu9bufD\nwB8BDuZN861XEamTEqucuHV/g7UPQ2oLOPOZ2kcBl5XBt7/tt3/4Q8jOrrmeRMcdd/iJIx57DObN\nq73e4J9BxllQshXe+5JewRGphxKrnJhP5/pRwAD590HHU2uv+6c/+XdXc3Phe9+LSnhSh7594fvf\n99vf+Ia/8alJShqM+Se0OAk+nQPzb6h5vmERAZRY5UTsXgJvTYKKo9DvJuhVx8o0H38Mt93mt//0\nJ98VLLF3223+Rmfp0srnrjVplQlnvwSprfzI7+UayCRSGyVWOT4H1sPcCZWDlUbcVftz1bIyvxzc\n4cP+9RrNCRw/WrWCv/7Vb//iF7BoUe11O4+EM5/2Ux6u+CV88pfoxCjSxCixSsPtXwtzxsHhbdD1\nPDjj7/4f29r89Kf+Gd7JJ/tRqBJfzj8fvvUtfwM0bRrs21d73cxLYGRgNaIFM+DjP0cnRpEmRIlV\nGmb3EvjPGDi4HjqNhLHPQWrz2uu/+KJ/TzI1Ff75T+jYMXqxSvjuuAMGDfLvFV9zjX8dqjZ9bqhc\nYm7hN/0yc3rmKvIZJVYJ3/a34fWz4fCn0HUcjJtT+whggA8/hC99yW/ffjucfXZ04pSGa9kSnnsO\nOnSA55/3vQx1GfC9wLSHBktvg8Xf1+xMIgFKrFI/52D1H333b+k+yL4cznkZmrWt/ZyCArj4Yjhw\nwHcvahRw/OvTx/cqpKbCr39d+wo4n9Wf7ieQsDRY/XuYe5HecxVBiVXqc3Q3vD0ZPvwuuDK/pNiY\np+ru/l29GsaOhU8/9c/vHnnET18o8W/8eHj4Yb/93e/W/0w890o477XKV3FeGQ7F7zV+nCJxTP/a\nSc2cg03PwMuDoWgmNGsPZz0HI/4AKam1nzdvnk+qW7fCOef4RbXT06MWtkTANdfAPff47e98x8+Q\nVdcz167nwvgPIWMMlBTB62Nh4bd874ZIElJilWPt+wTenODn/i0pgs6jYMJiyJ5c+znOwYMP+mRa\nXOxbPi+/DG3r6C6W+HXTTZU9Db/+NUyaVPtKOODfcx03F/JuBcwvdP9Sf9j4L79oukgSUWKVSntX\nw/vXwMsDYOtsaNYBRv4ZLpgHbepY2m3rVpg8Gb72NSgt9a9uvPiif0dSmq6vfMXfHHXsCC+9BIMH\nwyuv1F4/pRkM/Q2MXwSdTwtMgTgVXhnqez+UYCVJKLEmu4py2PIKvDMFXs6D9X/3+3tdB59bA32/\nXnvX76FD8NvfQr9+MHMmtG8Pjz/un8ulpUXvd5DGM368nzTitNOgqMgPSPv85/1rObXpOAQunAej\nHoSWmbBnue/9eHkQrL4LjuyMXvwiMaD1WJNRRTns/AAKZ8KGJ3x3L/gWR6/rIO8H0KaO5dy2bfMD\nXO6+23f7Akyc6KcqzMlp/Pgl+srL/Sjhn/8cDh70XcSXX+6fwZ5+eu2zbpUfgXWPwMrb4VCh35eS\nDpkTIfvz0H0CpOvd5kjQeqzxI6zEambjgbuBVOAh59wd1Y5b4PjFwCHgWufch+GcWxMl1girKPOt\nhh3z/IjNba9VbTW06e3n+e11rX9WVpNt22D2bL/k2+zZlRO2jxwJv/wlXHRRo/8aEge2boWf/cyv\niFMaWOVmwAC44gr/HHbIkJpHgFeUQtGLUPAQbH0VCPy7Y6l+0FPGWMg4E7qcDunto/XbJBQl1vhR\nb2I1s1TgY+ACoBBYAFzpnPsopM7FwE34xHoacLdz7rRwzq2JEutxcM4ny4Mb/KxIB9bD3o9g7wpf\nlpdUrd+mN2ReCtmX+SXBgi0O52DnTt/Vt3y57wZ8/31YubLy3JQU30L9xjf86zS1tVYkcRUWwr33\n+gT76aeV+zt3hjFj/HqvQ4ZAXp7vxQh9NHBwM2x+1o823/4WuGoTS7Tu6VdJaj8Q2vTxvSdtcqHF\nyZCqEea1UWKNH+Ek1jOAnzvnLgp8/iGAc+43IXX+ArzpnPtn4PMa4Bwgt75za3JcibWiFI7uCaPi\ncXR91/pn5Cqv50K2cSHnuMCgjUDpKoBA6corfypK/XuiFaW++6ziKFQc8QmxvATKSqDsgJ/0vnS/\nf7/06C6fTA9/6uftrWmdzHXAUSDtZGhxCjTvA816QmlrPyfs7t2wY4f/x3HLFti0qea5Ylu29DMn\nTZrkn7F17drwP0dJPKWlMGeOn7XplVdg8+Zj66Sl+bV3s7KgWzfIyPADojp2hFapULYJytZD542Q\nvtL/3a9N887Qohs07wLpnaB5J0hr52cAa9YW0lpBaku/Ck9qc9/tnNLcL31nzfzjDkv1PylpQErg\nc0rgBjFYBn4+2ybkM5Wfa9o+rhvNkHPSOwZia+AVlFjjRjj/9TKB0P9bCvGt0vrqZIZ5bmTsXgKz\nRzXKpZuMZu2hdY6/42/TE9r1h6/9Dxw4BGwN/LxZ/3XatvVrdQ4aBEOH+mdow4dD8zomhZDk1KyZ\nH+A0fry/mVy7FubPhwULYMUKP1lIURGsX+9/6vLAA3DVdbD/E9i9DPZ95HteDqzzvTCHP/U3kok+\n+OmSVdC+f6yjkBMQN0M3zWw6MB0g53gGwFiav4sNr/JxXL+2c0LuaGu9u00JrP5iVe+OLS3k7rmZ\n/5ySFrjDTvd33KktK3/S2gTuytv6u9rmnfxde4uu/iethjVOz3jWj95t3hxatPAtzzZt/E/79r7V\n0Lmzb0l06+a77Tp1UveuNJyZnxaxTx8/jWVQSYnvOi4shO3b/c+uXbB3r+8dOXjQ//Tu7f8/aJ/n\nf6qrKIcjxb535sgu32NzdJfvwSnd58vyEig/FCiPBHp9jlT2BlWUhvQUlVX2IrlyPuuB+qyHqbZe\nKEL2V9s+rsGg1c6xOiZgkSYhnMRaBGSHfM4K7AunTrMwzgXAOfcg8CD4ruAw4qqq0zCYUtzg0xLe\na6/FOgJJdi1b+h6Qvn1P7DopqdCym/8RiWPhvMe6AOhrZj3NLB2YCsysVmcmcLV5pwN7nXNbwzxX\nREQkYdTbYnXOlZnZjcBs/CszjzjnVprZjMDxB4BZ+BHBBfjXbb5S17mN8puIiIjEAU0QISKSADQq\nOH5oSkMREZEIUmIVERGJICVWERGRCFJiFRERiSAlVhERkQiKy1HBZlYMbIx1HAFdgB2xDqKBFHN0\nKOboaYpxRzvmHs65jCh+n9QiLhNrPDGzhU1tCLtijg7FHD1NMe6mGLNEhrqCRUREIkiJVUREJIKU\nWOv3YKwDOA6KOToUc/Q0xbibYswSAXrGKiIiEkFqsYqIiESQEmsYzOwmM1ttZivN7P/FOp5wmdn/\nmJkzs3BXgI8pM/td4M95mZk9Z2YdYh1TbcxsvJmtMbMCM7s11vHUx8yyzWyumX0U+Hv87VjHFC4z\nSzWzxWb2UqxjCYeZdTCzZwJ/l1eZ2RmxjkmiS4m1HmZ2LjAJGOKcGwjcGeOQwmJm2cCFwKZYx9IA\n/wEGOedOBT4GfhjjeGpkZqnAfcAEIA+40szyYhtVvcqA/3HO5QGnA99sAjEHfRtYFesgGuBu4FXn\nXH9gCE0rdokAJdb6fR24wzl3BMA5tz3G8YTrj8APgCbzEN0595pzrizw8QMgK5bx1GEUUOCcW+ec\nOwo8hb/5ilvOua3OuQ8D2/vx/9hnxjaq+plZFnAJ8FCsYwmHmbUHxgIPAzjnjjrn9sQ2Kok2Jdb6\n9QPOMrP5ZvaWmY2MdUD1MbNJQJFzbmmsYzkB1wGvxDqIWmQCm0M+F9IEklSQmeUCw4D5sY0kLHfh\nbxArYh1ImHoCxcCjge7rh8ysdayDkuhKi3UA8cDMXge61XDoR/g/o0747rORwNNm1svFeDh1PTHf\nhu8Gjjt1xe2ceyFQ50f4rssnohlbMjCzNsC/ge845/bFOp66mNmlwHbn3CIzOyfW8YQpDRgO3OSc\nm29mdwO3Aj+JbVgSTUqsgHPu/NqOmdnXgWcDifS/ZlaBnwO0OFrx1aS2mM1sMP6ueamZge9O/dDM\nRjnntkUxxBrV9WcNYGbXApcC42J981KHIiA75HNWYF9cM7Nm+KT6hHPu2VjHE4YxwEQzuxhoAbQz\ns8edc1fFOK66FAKFzrlgb8Az+MQqSURdwfV7HjgXwMz6AenE8WTgzrnlzrmTnHO5zrlc/P/ow+Mh\nqdbHzMbju/0mOucOxTqeOiwA+ppZTzNLB6YCM2McU53M32U9DKxyzv0h1vGEwzn3Q+dcVuDv8VTg\njThPqgT+P9tsZqcEdo0DPophSBIDarHW7xHgETNbARwFronjllRTdy/QHPhPoLX9gXNuRmxDOpZz\nrszMbgRmA6nAI865lTEOqz5jgC8Dy81sSWDfbc65WTGMKVHdBDwRuOlaB3wlxvFIlGnmJRERkQhS\nV7CIiEgEKbGKiIhEkBKriIhIBCmxioiIRJASq4iISAQpsYqIiESQEquIiEgEKbGKiIhE0P8H0FwT\ncnQiYkUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11500c18>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mu, sigma = 0, 1\n",
    "mu3, sigma3 = 0, (1 + 24*0.026)**0.5\n",
    "s = np.random.normal(mu, sigma, 10000)\n",
    "s3 = np.random.normal(mu, sigma3, 10000)\n",
    "#plot3 = count3, bins3, ignored3 = plt.hist(s3, 100, normed=True, color='y')\n",
    "#plot1 = count, bins, ignored = plt.hist(s, 100, normed=True, color='b')\n",
    "plt.plot(bins3, 1/(sigma3 * np.sqrt(2 * np.pi)) * np.exp( - (bins3 - mu)**2 / (2 * sigma3**2) ), linewidth=2, color='orange', label='mixture')\n",
    "plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) * np.exp( - (bins - mu)**2 / (2 * sigma**2) ), linewidth=2, color='r', label='normal')\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As demonstrated with the lines of fit for each Gaussian distribution, the Gaussian mixture has much thicker tails and a lower amplited. This would make sense, as the Gaussian distribution with a variance of 25 has far thicker tails and would spread the distribution among more points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1\n",
    "Find the Fisher information If associated with the Cauchy distribution, C(0,a), which has a density function given by\n",
    "$$ f(r) = \\frac{1}{\\pi a} \\frac{1}{1+(r/a)^{2}} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computation:\n",
    "$$ \\rho(r) = -ln(\\pi a) - ln(1 + (r/a)^{2}) $$\n",
    "Derive $\\rho(r)$ to find $\\psi(r) $:\n",
    "$$ \\psi(r) = \\frac{2(r/a)}{1+(r/a)^{2}}$$\n",
    "Derive $\\psi(r)$ to find $\\psi'(\\theta) $:\n",
    "$$ \\psi'(r) = -2\\frac{1-(r/a)^{2}}{(1+(r/a)^{2})^{2}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While Fisher information has been defined in our class notes as $I_{f}(\\theta) = E\\{[\\partial(lnf(z:\\theta))]^{2}\\}$ a preferred form, which is equivalent if the partial derivative exists, is actually $I_{f}(\\theta) = E\\{-[\\partial^{2}(lnf(z:\\theta))]\\}$ so:\n",
    "$$ I(r) = E\\{-[\\partial^{2}(lnf(r))]\\} = E\\{-[\\psi'(r)]\\} $$\n",
    "\n",
    "$$ = 2* \\int^{\\infty}_{-\\infty} 2\\frac{1-(r/a)^{2}}{(1+(r/a)^{2})^{2}} * f(r)dr $$\n",
    "\n",
    "$$ = \\frac{1}{\\pi a} \\int^{\\infty}_{-\\infty} 2\\frac{1-(r/a)^{2}}{(1+(r/a)^{2})^{2}}\\frac{1}{1+(r/a)^{2}} dr$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ = \\frac{2}{\\pi a} \\int^{\\infty}_{-\\infty} \\frac{1-(r/a)^{2}}{(1+(r/a)^{2})^{3}} dr$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Entering this integral into computation while replacing a with any real integer that is greater than 0 will result in:\n",
    "$$ I_{f}(r/a) = \\frac{1}{2} \\text{ for all real values of a > 0}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1\n",
    "At the Cauchy distribution C(0,1) with a = 1, calculate the asymptotic variance of the sample median, which is expressed as:\n",
    "$$ lim_{m\\rightarrow\\infty}V(\\sqrt{m}\\theta_{m}) = \\frac{E[\\psi^{2}(r)]}{(E[\\psi'(r)])^{2}} $$\n",
    "where E[.] is evaluated at the probability density function f(r) given with a = 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ f(r) = \\frac{1}{\\pi} \\frac{1}{1+(r)^{2}} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computation to help later:\n",
    "$$ \\rho(r) = -ln(\\pi) - ln(1 + (r)^{2}) $$\n",
    "Derive $\\rho(r)$ to find $\\psi(r) $:\n",
    "$$ \\psi(r) = \\frac{2(r)}{1+(r)^{2}}$$\n",
    "Square the $ \\psi(r) $ function for later use:\n",
    "$$ \\psi^{2}(r) = (\\frac{2(r)}{1+(r)^{2}})^{2} = \\frac{4r^{2}}{(1+r^{2})^{2}}$$\n",
    "Derive $\\psi(r)$ to find $\\psi'(r) $:\n",
    "$$ \\psi'(r) = -2\\frac{1-(r)^{2}}{((1+(r)^{2})^{2}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since $I(r) = E\\{-[\\partial^{2}(lnf(r))]\\} = E\\{-[\\psi'(r)]\\} = 1/2$ (derived in the previous problem) then $[E\\{-[\\psi'(r)]\\}]^{2} = E\\{[\\psi'(r)]\\}^2 = 1/4$\n",
    "\n",
    "$$ E[\\psi^{2}(r)] = \\int^{\\infty}_{-\\infty} \\frac{4r^{2}}{(1+r^{2})^{2}} * \\frac{1}{\\pi} \\frac{1}{1+(r)^{2}} dr = \\frac{4}{\\pi}\\int^{\\infty}_{-\\infty} \\frac{r^{2}}{(1+r^{2})^{3}} dr$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Entering the integral above into computation results in a value of 1/2 so\n",
    "$$ lim_{m\\rightarrow\\infty}V(\\sqrt{m}\\theta_{m}) = \\frac{E[\\psi^{2}(r)]}{(E[\\psi'(r)])^{2}} = \\frac{1/2}{1/4} = 2 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2\n",
    "Calculate the asymptotic efficiency of the sample median and its loss of efficiency:\n",
    "\n",
    "Asymptotic efficiency is defined by: \n",
    "$$e_{m} = \\frac{1/I_{f}(r)}{V(\\sqrt{m}\\theta_{m})} = \\frac{\\frac{1}{1/2}}{2} = 2/2 = 1$$\n",
    "\n",
    "This would indicate that the variance of the sample mean is an **efficient estimator** as it asymptotically approaches $m\\rightarrow \\infty$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1\n",
    "At the standard Gaussian distribution, N(0,1), calculate the asymptotic variance given by (1) of the Huber estimator (b=1.5) and the Huber-type skipped mean (b=2.71).\n",
    "$$ lim_{m\\rightarrow\\infty}V(\\sqrt{m}\\theta_{m}) = \\frac{E[\\psi^{2}(r)]}{(E[\\psi'(r)])^{2}} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For both estimators:\n",
    "\n",
    "$f(r) = \\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}} \\ \\ for \\ \\ N(0,1) $\n",
    "\n",
    "For the Huber-Estimator:\n",
    "\n",
    "$\\begin{equation}\n",
    "    \\rho(r)=\n",
    "    \\begin{cases}\n",
    "      r^{2}/2, & \\text{for}\\ |r|\\leq b \\\\\n",
    "      b|r| - \\frac{b^{2}}{2}, & \\text{for}\\ |r|>b \\\\\n",
    "    \\end{cases}\n",
    "  \\end{equation}$\n",
    "\n",
    "$\\begin{equation}\n",
    "\\psi(r)=\n",
    "\\begin{cases}\n",
    "  r, & \\text{for}\\ |r|\\leq b \\\\\n",
    "  bsign(r), & \\text{for}\\ |r|>b \\\\\n",
    "\\end{cases}\n",
    "\\end{equation}$\n",
    "\n",
    "For the Huber-type skipped mean:\n",
    "  \n",
    "$\\begin{equation}\n",
    "\\rho(r)=\n",
    "\\begin{cases}\n",
    "  r^{2}/2, & \\text{for}\\ |r|\\leq b \\\\\n",
    "  \\frac{b^{2}}{2}, & \\text{for}\\ |r|>b \\\\\n",
    "\\end{cases}\n",
    "\\end{equation}$\n",
    "\n",
    "$\\begin{equation}\n",
    "\\psi(r)=\n",
    "\\begin{cases}\n",
    "  r, & \\text{for}\\ |r|\\leq b \\\\\n",
    "  0, & \\text{for}\\ |r|>b \\\\\n",
    "\\end{cases}\n",
    "\\end{equation}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we've stated above:\n",
    "$$ lim_{m\\rightarrow\\infty}V(\\sqrt{m}\\theta_{m}) = \\frac{E[\\psi^{2}(r)]}{(E[\\psi'(r)])^{2}} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**For Huber-Estimator:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ E[\\psi^{2}(2)] = \\int^{\\infty}_{-\\infty}\\psi^{2}(r)f(r)dr = $$\n",
    "$$\\int^{-1.5}_{-\\infty}(bsign(r))^{2}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr + \\int^{1.5}_{-1.5}r^{2}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr + \\int^{\\infty}_{1.5}(bsign(r))^{2}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since $(sign(x)^2) = 1$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$(b)^{2}\\int^{-1.5}_{-\\infty}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr + \\int^{1.5}_{-1.5}r^{2}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr + (b)^{2}\\int^{\\infty}_{1.5}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr$$\n",
    "$$ E[\\psi^{2}(2)] = (1.5)^{2}\\int^{-1.5}_{-\\infty}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr + \\int^{1.5}_{-1.5}r^{2}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr + (1.5)^{2}\\int^{\\infty}_{1.5}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computing the full integration, you get: $E[\\psi^{2}(2)] = 0.150316 + 0.477833 + 0.150316 = 0.778465 $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\begin{equation}\n",
    "\\psi'(r)= \\partial \\psi(r) = \n",
    "\\begin{cases}\n",
    "  1, & \\text{for}\\ |r|\\leq b \\\\\n",
    "  b2\\delta(r), & \\text{for}\\ |r|>b \\\\\n",
    "\\end{cases}\n",
    "\\end{equation} \\ \\ \\ $\n",
    "$Where \\ \\delta(r) = Dirac Function = \n",
    "\\begin{equation}\n",
    "\\begin{cases}\n",
    "  +\\infty, & \\text{for}\\ r = 0 \\\\\n",
    "  0, & \\text{for}\\ r \\neq 0 \\\\\n",
    "\\end{cases}\n",
    "\\end{equation}$ so effectively 0 for this problem.\n",
    "$$ E[\\psi'(r)] = $\\int^{\\infty}_{-\\infty}\\psi'(r)f(r)dr = $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\int^{-1.5}_{-\\infty}0^{2}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr + \\int^{1.5}_{-1.5}1^{2}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr + \\int^{\\infty}_{1.5}0^{2}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr = \\int^{1.5}_{-1.5}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computing the full integration you get: $E[\\psi'(r)] = 0.866386$; $E[\\psi'(r)]^2 = 0.750624$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ lim_{m\\rightarrow\\infty}V(\\sqrt{m}\\theta_{m}) = \\frac{E[\\psi^{2}(r)]}{(E[\\psi'(r)])^{2}} = \\frac{0.778465}{0.750624} = 1.03709 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**For Huber-type skipped mean:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\begin{equation}\n",
    "\\psi(r)=\n",
    "\\begin{cases}\n",
    "  r, & \\text{for}\\ |r|\\leq b \\\\\n",
    "  0, & \\text{for}\\ |r|>b \\\\\n",
    "\\end{cases}\n",
    "\\end{equation}$\n",
    "\n",
    "$\\begin{equation}\n",
    "\\psi'(r)= \\partial \\psi(r) = \n",
    "\\begin{cases}\n",
    "  1, & \\text{for}\\ |r|\\leq b \\\\\n",
    "  0, & \\text{for}\\ |r|>b \\\\\n",
    "\\end{cases}\n",
    "\\end{equation} \\ \\ \\ $\n",
    "$$ E[\\psi^{2}(2)] = \\int^{\\infty}_{-\\infty}\\psi^{2}(r)f(r)dr = $$\n",
    "$$\\int^{-1.5}_{-\\infty}(0)^{2}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr + \\int^{1.5}_{-1.5}r^{2}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr + \\int^{\\infty}_{1.5}(0)^{2}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computing the full integration you get $E[\\psi^{2}(2)] = 0.477833$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$E[\\psi'(r)] = \\int^{-1.5}_{-\\infty}0^{2}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr + \\int^{1.5}_{-1.5}1^{2}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr + \\int^{\\infty}_{1.5}0^{2}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr = \\int^{1.5}_{-1.5}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computing the full integration you get: $E[\\psi'(r)] = 0.866386$; $E[\\psi'(r)]^2 = 0.750624$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ lim_{m\\rightarrow\\infty}V(\\sqrt{m}\\theta_{m}) = \\frac{E[\\psi^{2}(r)]}{(E[\\psi'(r)])^{2}} = \\frac{0.477833}{0.750624} = 0.636581$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2\n",
    "calculate the asymptotic relative efficiency and the loss of efficiency of each of these two estimators with respect to the sample mean."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**For sample mean of Gaussian N(0,1)**\n",
    "$$ I_{f}(r) = E[\\psi^{2}(r)] = \\int^{\\infty}_{-\\infty} r^{2}f(r)dr = \\int^{\\infty}_{-\\infty}r^{2}\\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{(r)^{2}}{2}}dr =1$$\n",
    "$$ lim_{m\\rightarrow \\infty} V(\\sqrt{m}\\bar{\\theta}_{m}) = \\frac{1}{I_{f}} = 1$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**For Huber Estimator:**\n",
    "Asymptotic relative efficiency: \n",
    "$$e_{\\theta_{\\mu}\\theta{he}} = lim_{m\\rightarrow\\infty}\\frac{V(\\sqrt{m}\\theta_{\\mu})}{{V(\\sqrt{m}\\theta_{\\mu})}} = \\frac{1}{1.03709} = 0.964236 $$\n",
    "There is a loss of efficiency of 3.4%, making the Huber estimator with b = 1.5 a very good trade-off for efficiency in exchange for robust-ness of including outliers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**For Huber Skipped-mean:**\n",
    "Asymptotic relative efficiency: \n",
    "$$e_{\\theta_{\\mu}\\theta{he}} = lim_{m\\rightarrow\\infty}\\frac{V(\\sqrt{m}\\theta_{\\mu})}{{V(\\sqrt{m}\\theta_{\\mu})}} = \\frac{1}{0.636581} = 1.57 \\approx \\pi/2$$\n",
    "There is a loss of efficiency of 57%, making the Huber estimator with b = 1.5 incredibly impractical and exchanges too much efficiency for any amount of robust-ness it may add, at least without contamination of the sample."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}