{ "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, 2nd Edition](https://www.crcpress.com/Introduction-to-Probability-Second-Edition/Blitzstein-Hwang/p/book/9781138369917), 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 6th 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−9. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "☣ 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 nth 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 2nd 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 2nd 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 nth 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 6th sample moment of 100 i.i.d. $N(0, 1)$ r.v.s. How close is it to the true 6th 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", "☣ 6.9.2. scipy.optimize.rootis 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/66. 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/66 $\\approx$ 0.07354." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "\n", "Joseph K. Blitzstein and Jessica Hwang, Harvard University and Stanford University, © 2019 by Taylor and Francis Group, LLC" ] } ], "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 }