Ch6.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 6: Moments\n",
    " \n",
    "This Jupyter notebook is the Python equivalent of the R code in section 6.9 R, pp. 267 - 272, [Introduction to Probability, 1st Edition](https://www.crcpress.com/Introduction-to-Probability/Blitzstein-Hwang/p/book/9781466575578), Blitzstein & Hwang.\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions\n",
    "\n",
    "The MGF of an r.v. is a _function_. As an example of defining and working with functions in Python, let's use the $N(0, 1)$ MGF, which is given by $M(t) = e^{\\frac{t^{2}}{2}}$. The code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calling function M with a single value: \n",
      "M(0) = 1.0\n",
      "\n",
      "calling M with a vector: \n",
      "M(np.arange(1,11))) = [1.64872127e+00 7.38905610e+00 9.00171313e+01 2.98095799e+03\n",
      " 2.68337287e+05 6.56599691e+07 4.36731791e+10 7.89629602e+13\n",
      " 3.88084696e+17 5.18470553e+21]\n"
     ]
    }
   ],
   "source": [
    "def M(t):\n",
    "    \"\"\" Return N(0, 1) MGF evalated at t. \"\"\"\n",
    "    return np.exp(t**2/2)\n",
    "\n",
    "print('calling function M with a single value: \\nM(0) = {}\\n'.format(M(0)))\n",
    "\n",
    "print('calling M with a vector: \\nM(np.arange(1,11))) = {}'.format(M(np.arange(1,11))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "defines `M` to be this function. The `def` Python keyword says that we're [defining a function](https://docs.python.org/3/tutorial/controlflow.html#defining-functions). The function is named `M`, and it takes one variable `t` (called the argument or parameter of the function). The line declaring the function name and list of parameter(s) is terminated with a colon `:`, with the body of the function following on the next line after an indent. Note that a simple Python function will not be able to flexibly deal with both single or vector inputs, such as is possible with a function in R. However, since our function `M` relies on [`numpy.exp`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.exp.html), `M` can accept both a nested sequence of objects or `numpy` arrays as inputs, and return an single or tuple of `numpy` array as output."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Writing\n",
    "\n",
    "    def M(x):\n",
    "        return np.exp(x**2/2)\n",
    "\n",
    "would define the same function `M`, except that now the parameter is named `x`. Giving the parameters names is helpful for functions of more than one parameter, since Python then saves us from having to remember the order in which to write the parameters, and allows us to assign default values. For example, the $N(\\mu, \\sigma^2)$ MGF is given by $g(t) = exp\\left(\\mu \\, t + \\frac{\\sigma^2 \\, t^2}{2} \\right)$, which depends on $t$, $mu$, and $\\sigma$. We can define this in Python by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N(2, 3) MGF evaluated at 1 = 665.1416330443618\n"
     ]
    }
   ],
   "source": [
    "def g(t, mean=0, sd=1):\n",
    "    \"\"\" Return the N(mean, sd) MGF evaluated at t.\n",
    "        default mean (mu) = 0\n",
    "        default sd (sigma) = 1\n",
    "    \"\"\"\n",
    "    return np.exp(mean*t + (sd**2 * t**2)/2)\n",
    "\n",
    "ans = g(1, 2, 3)\n",
    "print('N(2, 3) MGF evaluated at 1 = {}'.format(ans))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is `g(1, 2, 3)`? It's the $N(2, 3^{2})$ MGF evaluated at 1, but it may be hard to remember which parameter is which, especially when working with many functions with many parameters over the course of many months. So we can also write `g(1, mean=2, sd=3)` or `g(1, sd=3, mean=2)`. Since the `mean` and `sd` function parameters have the form `parameter = expression`, the function is said to have \"default parameter values.\"\n",
    "\n",
    "Also, when defining `g` we specified default values of 0 for the `mean` and 1 for the `sd` standard deviation, so if we want the $N(0, 5^2)$ MGF evaluated at 3, we can use `g(3, sd=5)` as shorthand. It would be bad here to write `g(3, 5)`, since that is ambiguous about which argument is omitted; in fact, Python interprets this as `g(t3, mean=5)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "g(1, mean=2, sd=3) = 665.1416330443618\t... explicitly using parameter names\n",
      "g(3, 5) = 294267566.0415088\t\t... but which parameter was omitted?\n",
      "g(3, mean=5) = 294267566.0415088\t... 'mean' parameter was omitted\n"
     ]
    }
   ],
   "source": [
    "ans1 = g(1, mean=2, sd=3)\n",
    "print('g(1, mean=2, sd=3) = {}\\t... explicitly using parameter names'.format(ans1))\n",
    "\n",
    "ans2 = g(3, 5)\n",
    "print('g(3, 5) = {}\\t\\t... but which parameter was omitted?'.format(ans2))\n",
    "\n",
    "ans3 = g(3, mean=5)\n",
    "print('g(3, mean=5) = {}\\t... \\'mean\\' parameter was omitted'.format(ans3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Moments\n",
    "\n",
    "LOTUS makes it easy to write down any moment of a continuous r.v. as an integral. The [`scipy.integrate`](https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html) module in SciPy can help us do the integral numerically, using the [`scipy.integrate.quad`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html#scipy.integrate.quad) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import quad\n",
    "\n",
    "# to learn more about scipy.integrate.quad, un-comment out the following line\n",
    "#print(quad.__doc__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, let's approximate the 6<sup>th</sup> moment of a $N(0, 1)$ r.v. The code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6th moment of N(0,1) = 15.000000000000004, with error = 4.4229319318339374e-09\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "def g(x):\n",
    "    \"\"\" Return the 6th moment of N(0, 1). \"\"\"\n",
    "    return x**6 * norm.pdf(x)\n",
    "\n",
    "y, abserr = quad(g, -np.inf, np.inf)\n",
    "print('6th moment of N(0,1) = {}, with error = {}'.format(y, abserr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "asks `quad` to compute $\\int_{-\\infty}^{\\infty} g(x) \\, dx$, where $g(x) = x^6 \\, \\phi(x)$ with $\\phi$ the $N(0, 1)$ PDF. When we ran this, `quad` reported 15 (the correct answer, as we know from this chapter!) and that the absolute error was less than 4.423 $\\times$ 10<sup>−9</sup>. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&#x2623; 6.9.1. Numerical integration runs into difficulties for some functions; as usual, checking answers in multiple ways is a good idea. Using `numpy.inf` for parameter `b` (the upper limit of integration) is preferred to using a large number as the upper limit when integrating up to $\\infty$ (and likewise for a lower limit of `a = -numpy.inf` for $-\\infty$. For example, on many systems `quad(norm.pdf, 0, 10**6)` reports 0.0 while `quad(norm.pdf, 0, numpy.inf)` reports the correct answer, 0.5."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, all of the continuous distributions supported in [`scipy.stats`](https://docs.scipy.org/doc/scipy/reference/tutorial/stats/continuous.html#moments) have a `moment(n, loc=0, scale=1)` function that allows you to quickly and easily obtain the n<sup>th</sup> order non-central moment of the continuous distribution in question. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "norm.moment(n=6) = 15.00000000089533\n"
     ]
    }
   ],
   "source": [
    "print('norm.moment(n=6) = {}'.format(norm.moment(n=6)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, to check that the 2<sup>nd</sup> moment (and variance) of a $Unif(−1, 1)$ r.v. is 1/3, we can use `quad` again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2nd moment of Unif(-1,1) = 0.3333333333333333, with error = 3.700743415417188e-15\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import uniform\n",
    "\n",
    "def h(x):\n",
    "    \"\"\" Return 2nd moment (var) of Unif(-1, 1).\n",
    "    \n",
    "        scipy.stats.uniform is constant between\n",
    "            a = loc\n",
    "            b = loc + scale\n",
    "    \"\"\"\n",
    "    a = -1\n",
    "    b = 1\n",
    "    loc = a\n",
    "    scale = b - loc\n",
    "    return x**2 * uniform.pdf(x, loc=loc, scale=scale)\n",
    "\n",
    "y, abserr = quad(h, -1, 1)\n",
    "print('2nd moment of Unif(-1,1) = {}, with error = {}'.format(y, abserr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternately, we can either use `uniform.moment(n=2, loc=-1, scale=2)` or just `uniform.var(loc=-1, scale=2)`, keeping in mind that `loc=-1` and `scale = interval length = 2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "uniform.moment(n=2, loc=-1, scale=2) = 0.33333333333333326\n",
      "uniform.var(loc=-1, scale=2) = 0.3333333333333333\n"
     ]
    }
   ],
   "source": [
    "ans1 = uniform.moment(n=2, loc=-1, scale=2)\n",
    "print('uniform.moment(n=2, loc=-1, scale=2) = {}'.format(ans1))\n",
    "\n",
    "ans2 = uniform.var(loc=-1, scale=2)\n",
    "print('uniform.var(loc=-1, scale=2) = {}'.format(ans2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For moments of a discrete r.v., we can use LOTUS and the `numpy.sum` function. For example, to find the 2<sup>nd</sup> moment of $X \\sim Pois(7)$, we can use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2nd moment of Pois(7) = 55.99999999999994\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import poisson\n",
    "\n",
    "def g(k):\n",
    "    \"\"\" Return the 2nd moment of Pois(7) evaluated at k. \"\"\"\n",
    "    return k**2 * poisson.pmf(k, 7)\n",
    "\n",
    "# we want to sum up to and including 100, so the upper limit is 100+1\n",
    "ans = np.sum(g(np.arange(0, 100+1)))\n",
    "print('2nd moment of Pois(7) = {}'.format(ans))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we summed up to 100 since it’s clear after getting a sense of the terms that the total contribution of all the terms after k = 100 is negligible (choosing an upper limit in this way is in contrast to how to use the integrate command in the continuous case). The result is extremely close to 56, which is comforting since $E(X^{2}) = Var(X) + (EX)^{2} = 7 + 49 = 56$.\n",
    "\n",
    "But similar to continous r.v., the [discrete r.v. in `scipy.stats`](https://docs.scipy.org/doc/scipy/reference/tutorial/stats/discrete.html#moments) have a `moment` function as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "poisson.moment(n=2, mu=7) = 56\n"
     ]
    }
   ],
   "source": [
    "ans = poisson.moment(n=2, mu=7)\n",
    "print('poisson.moment(n=2, mu=7) = {}'.format(ans))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A sample moment can be found in easily using NumPy. If `x` is a vector of data, then [`numpy.mean(x)`](https://docs.scipy.org/doc/numpy-1.15.1/reference/generated/numpy.mean.html) gives its sample mean and, more generally, `numpy.mean(x**n)` gives the n<sup>th</sup> sample mean for any positive integer `n`. For example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sample moment: numpy.mean(x**6) = 9.336344243424868\n"
     ]
    }
   ],
   "source": [
    "# seed the random number generator\n",
    "np.random.seed(6765)\n",
    "\n",
    "x = norm.rvs(size=100)\n",
    "print('sample moment: numpy.mean(x**6) = {}'.format(np.mean(x**6)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "gives the 6<sup>th</sup> sample moment of 100 i.i.d.  $N(0, 1)$ r.v.s. How close is it to the true 6<sup>th</sup> moment? How close are other sample moments to the corresponding true moments?\n",
    "\n",
    "The sample variance can also be found in easily with NumPy. If `x` is a vector of data, then using the `ddof` parameter (delta degrees of freedom) such as in [`numpy.var(x, ddof=1)`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.var.html) gives its sample variance. This returns `nan` (not a number) as well as issuing `RuntimeWarning: Degrees of freedom <= 0 for slice` (the divisor used in the calculation is `len(x) - ddof`) if `x` has length 1, since the $n − 1$ in the denominator is 0 in this case. It makes sense not to return a numerical value in this case, not only because of the definition but also because it would be insane to try to estimate the variability of a population if we only have one observation!\n",
    "\n",
    "For a simple demonstration of using the sample mean and sample variance to estimate the true mean and true variance of a distribution, we generate 1000 times from a $N(0, 1)$ distribution and store the values in `z`. We then compute the sample mean and sample variance with `numpy.mean(z)` and `numpy.var(z, ddof=1)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sample mean of z: 0.00915384129813555\n",
      "sample variance of z: 1.073340109262066\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(10946)\n",
    "\n",
    "z = norm.rvs(size=1000)\n",
    "\n",
    "mu_z = np.mean(z)\n",
    "print('sample mean of z: {}'.format(mu_z))\n",
    "\n",
    "var_z = np.var(z, ddof=1)\n",
    "print('sample variance of z: {}'.format(var_z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We find that `numpy.mean(z)` is close to 0 and `numpy.var(z, ddof=1)` is close to 1. You can try this out for a $N(\\mu, \\sigma^2)$ distribution (or other distribution) of your choosing; just remember that `numpy.norm.rvs` takes $\\sigma$ as the `scale` parameter, and not $\\sigma^2$!\n",
    "\n",
    "The sample standard deviation of `x` can be found using [`numpy.std(x, ddof=1)`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.std.html). This gives the same result as `numpy.sqrt(numpy.var(x, ddof=1))`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "np.std(z, ddof=1) = 1.0360212880351765\n",
      "np.sqrt(np.var(z, ddof=1)) = 1.0360212880351765\n"
     ]
    }
   ],
   "source": [
    "sd1 = np.std(z, ddof=1)\n",
    "print('np.std(z, ddof=1) = {}'.format(sd1))\n",
    "\n",
    "sd2 = np.sqrt(np.var(z, ddof=1))\n",
    "print('np.sqrt(np.var(z, ddof=1)) = {}'.format(sd2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While the [`scipy.stats`](https://docs.scipy.org/doc/scipy/reference/stats.html) module does have functions for [`skew`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.skew.html) and [`kurtosis`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.kurtosis.html), both functions rely on the _population standard deviation_ (with $n$ rather than $n-1$ in the denominator).\n",
    "\n",
    "However, we can easily define our own functions for _sample skewness_ and _sample kurtosis_ by using [`scipy.stats.moment`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.moment.html) and [`numpy.std(z, ddof=1)`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.std.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import moment\n",
    "\n",
    "# to learn more about scipy.stats.moment, un-comment out the following line\n",
    "#print(moment.__doc__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sample skew of z = 0.02333931072753422\n"
     ]
    }
   ],
   "source": [
    "def skew(x, use_sample_sd=True):\n",
    "    \"\"\" Return the skew of x.\n",
    "    \n",
    "        Default is to use sample standard deviation on the denominator,\n",
    "        yielding sample skew. \n",
    "        \n",
    "        Specifying use_sample_sd=False is the same as using\n",
    "        scipy.stats.skew(x).\n",
    "    \"\"\"\n",
    "    ddof = 1 if use_sample_sd==True else 0\n",
    "    return moment(x, 3) / np.std(x, ddof=ddof)**3\n",
    "\n",
    "print('sample skew of z = {}'.format(skew(z)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sample kurtosis of z = 0.16247912720035718\n"
     ]
    }
   ],
   "source": [
    "def kurt(x, use_sample_sd=True):\n",
    "    \"\"\" Return the excess kurtosis of x.\n",
    "    \n",
    "        Default is to use sample standard deviation on the denominator,\n",
    "        yielding sample kurtosis. \n",
    "        \n",
    "        Specifying use_sample_sd=False is the same as using\n",
    "        scipy.stats.kurtosis(x).\n",
    "    \"\"\"\n",
    "    ddof = 1 if use_sample_sd==True else 0\n",
    "    return moment(x, 4) / np.std(x, ddof=ddof)**4 - 3.0\n",
    "\n",
    "print('sample kurtosis of z = {}'.format(kurt(z)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Medians and modes\n",
    "\n",
    "To find the median of a continuous r.v. with CDF $F$, we need to solve the equation $F(x) = 1/2$ for $x$, which is equivalent to finding the root (zero) of the function $g$ given by $g(x) = F(x) − 1/2$. This can be done using [`scipy.optimize.root`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.root.html#scipy.optimize.root) in SciPy. For example, let's find the median of the Expo(1) distribution. The code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "    fjac: array([[-1.]])\n",
       "     fun: array([-5.55111512e-17])\n",
       " message: 'The solution converged.'\n",
       "    nfev: 9\n",
       "     qtf: array([-1.56348268e-11])\n",
       "       r: array([0.5000001])\n",
       "  status: 1\n",
       " success: True\n",
       "       x: array([0.69314718])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import expon\n",
    "from scipy.optimize import root\n",
    "\n",
    "# to learn more about scipy.optimize.root, un-comment out the following line\n",
    "#print(root.__doc__)\n",
    "\n",
    "def g(x):\n",
    "    \"\"\" Assuming F(x) = Expo(1),\n",
    "        define a function g(x) = F(x) - 1/2\n",
    "    \"\"\"\n",
    "    return expon.pdf(x) - 1/2\n",
    "\n",
    "# set our inital guess at a root to be 0\n",
    "root(g, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "asks `root` to find a root of the desired function, with an initial guess of 0. This returns an answer very close to the true answer of $log(2) ≈ 0.693$. Of course, in this case we can solve $1 − e^{−x} = \\frac{1}{2}$ directly without having to use numerical methods.\n",
    "\n",
    "&#x2623; 6.9.2. `scipy.optimize.root`is useful but there is no guarantee that it will find a root (`success` may be `False`). When using `scipy.optmize.root` please pay attention to the `status`, `success` and `message` in the result returned.\n",
    "\n",
    "An easier way to find the median of the $Expo(1)$ is to use `scipy.stats.expon.median`. The `median` function of returns the median of the distribution, for `expon` as well as all other continuous distributions in `scipy.stats`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "expon.median() = 0.6931471805599453\n"
     ]
    }
   ],
   "source": [
    "ans = expon.median()\n",
    "print('expon.median() = {}'.format(ans))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For finding the mode of a continuous distribution, we can use the [`fminbound`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fminbound.html#scipy.optimize.fminbound) function in SciPy. For example, let's find the mode of the $Gamma(6, 1)$ distribution, which is an important distribution that we will introduce in the next chapter. Its PDF is proportional to $x^{5} \\, e^{−x}$. Using calculus, we can find that the mode is at $x = 5$. Using [`scipy.optimize.fminbound`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fminbound.html#scipy.optimize.fminbound), we can find that the mode is very close to $x = 5$ as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mode of Gamma(6, 1) is approximately 5.000000043355882\n"
     ]
    }
   ],
   "source": [
    "def h(x):\n",
    "    \"\"\" Since fminbound MINIMIZES the argument function,\n",
    "        we will instead pass in the negative of the Gamma PDF\n",
    "        in order to obtain the MAXIMUM.\n",
    "    \"\"\"\n",
    "    def f(x):\n",
    "        \"\"\" Gamma PDF is proportional to this function.\n",
    "        \"\"\"\n",
    "        return x**5 * np.e**-x\n",
    "    return -f(x)\n",
    "\n",
    "from scipy.optimize import fminbound\n",
    "\n",
    "lower_bound = 0.0\n",
    "upper_bound = 20.0\n",
    "\n",
    "ans = fminbound(h, lower_bound, upper_bound)\n",
    "print('mode of Gamma(6, 1) is approximately {}'.format(ans))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we had wanted to _minimize_ instead of _maximize_, we could have put simply passed in `f` to `fminbound` without composing the nested function `h`.\n",
    "\n",
    "Next, let's do a discrete example of median and mode. An interesting fact about the $Bin(n, p)$ distribution is that if the mean $np$ is an integer, then the median and mode are also $np$ (even if the distribution is very skewed). To check this fact about the median for the $Bin(50, 0.2)$ distribution, we can use the following code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "median of Bin(50, 0.2) is 10\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import binom\n",
    "\n",
    "bin_cdf_vector = binom.cdf(np.arange(0,50), 50, 0.2)\n",
    "median = np.argmax(bin_cdf_vector>=0.5)\n",
    "print('median of Bin(50, 0.2) is {}'.format(median))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [`numpy.argmax`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.argmax.html) function finds the location of the maximum of a vector, giving the index of the first occurrence of a maximum. Since `True` is encoded as 1 and `False` is encoded as 0, the first maximum in `binom.cdf(np.arange(0,50), 50, 0.2) >= 0.5` is at the first value for which the CDF is at least 0.5. The return value of `numpy.argmax` on `binom.cdf(np.arange(0,50), 50, 0.2)` is 10, showing that the median is at 10."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, `numpy.argmax` on `binom.pmf(np.arange(0,50), 50, 0.2)` returns 10, showing that the mode is at 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mode of Bin(50, 0.2) is 10\n"
     ]
    }
   ],
   "source": [
    "bin_pmf_vector = binom.pmf(np.arange(0,50), 50, 0.2)\n",
    "mode = np.argmax(bin_pmf_vector)\n",
    "print('mode of Bin(50, 0.2) is {}'.format(mode))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sample median of a vector x of data can be found using [`numpy.median(x)`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.median.html). [`scipy.stats.mode`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.mode.html) can be used to find the sample mode, but it will only return the smallest mode in case there are ties. Instead, we can compose our own `mode` function, using [`numpy.where`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.where.html) and [`numpy.bincount`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.bincount.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mode(s) in r.v. vector x: [8]\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(17711)\n",
    "\n",
    "bin_rv_vector = binom.rvs(50, 0.2, size=50)\n",
    "\n",
    "def mode(x):\n",
    "    m = np.bincount(x).max()\n",
    "    return np.where(np.bincount(bin_rv_vector)==m)[0]\n",
    "\n",
    "ans = mode(bin_rv_vector)\n",
    "print('mode(s) in r.v. vector x: {}'.format(ans))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dice simulation\n",
    "\n",
    "In the starred Section 6.7, we showed that in rolling 6 fair dice, the probability of a total of 18 is 3431/6<sup>6</sup>. But the proof was complicated. If we only need an approximate answer, simulation is a much easier approach. And we already know how to do it! Here is the code for a million repetitions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "probability of total 18 with 6 fair dice (simulated) is 0.073304\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(28657)\n",
    "\n",
    "repetitions = 10**6\n",
    "# counter to track number of 18s seen\n",
    "c = 0\n",
    "\n",
    "for _ in range(repetitions):\n",
    "    if np.random.choice(np.arange(1,7), 6).sum() == 18:\n",
    "        # increment the counter\n",
    "        c += 1\n",
    "\n",
    "ans = c / repetitions\n",
    "print('probability of total 18 with 6 fair dice (simulated) is {}'.format(ans))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our simulation this yielded 0.073304, which is very close to 3431/6<sup>6</sup> $\\approx$ 0.07354."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "\n",
    "&copy; Blitzstein, Joseph K.; Hwang, Jessica. Introduction to Probability (Chapman & Hall/CRC Texts in Statistical Science)."
   ]
  }
 ],
 "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.6.3"
  },
  "notebook_info": {
   "author": "Joseph K. Blitzstein, Jessica Hwang",
   "chapter": "6. Moments",
   "creator": "buruzaemon",
   "section": "6.9",
   "title": "Introduction to Probability, 1st Edition"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}