{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Week 4: Mathematical expectation\n",
"\n",
"\n",
" #### [Back to main page](https://petrosyan.page/fall2020math3215)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Mean value of the hypergeometric distribution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The mean value of the hypergeometric distribution with parameters $(N,K,n)$ is equal to \n",
"\n",
"$$\\mu=\\frac{nK}{N}.$$\n",
"\n",
"As we discussed in class, the mean value was in a sense the center of the histogram: the point around which most of the histogram was concentrated. This is very vividly observed for the hypergeometric distribution. The plot bellow shows the mean value in blue for different values of $(N,K,n)$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "b8355df531f74e289ec832d171d24e3f",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"interactive(children=(IntSlider(value=80, description='N', min=1, readout_format=''), IntSlider(value=40, desc…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# nbi:hide_in\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from scipy.special import comb\n",
"from ipywidgets import interact, IntSlider\n",
"\n",
"def hypergeometric_pmf(N=80,K=40,n=30):\n",
" range_x = np.arange(max(0, n-(N-K)), min(n, K)+1, 1)\n",
"\n",
" def hyper_pmf(N,K,n,i):\n",
" pmf_val = comb(K, i, exact=True) * comb(N-K, n-i, exact=True) / comb(N, n, exact=True)\n",
" return pmf_val\n",
" mean = n * K / N\n",
"\n",
" pmf_values = np.array([hyper_pmf(N,K,n,i) for i in range_x])\n",
"\n",
" # plot setup\n",
" plt.figure(figsize=(14,7)) \n",
" plt.axhline(y=0, color='k')\n",
" plt.ylim(-0.02, max(np.max(pmf_values)+0.05, 0.2))\n",
" plt.xlim(-2,N+2)\n",
" plt.xticks(np.arange(0, N+1, 5))\n",
"\n",
" # PLotting with plt.bar instead of plt.hist works better when f(x) are knowwn\n",
" plt.scatter(np.array([mean]),np.zeros(1), color =\"blue\", s=200)\n",
" plt.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20)\n",
" plt.bar(range_x, pmf_values, width=1, color=(0.1, 0.1, 1, 0.1), edgecolor='blue', linewidth=1.3)\n",
" plt.title(\"Hypergeometric distribution\")\n",
" plt.figtext(0.75,0.8, \" N={} \\n K={} \\n n={} \\n expectation={:.2f}\".format(N,K,n, mean), ha=\"left\", va=\"top\",\n",
" backgroundcolor=(0.1, 0.1, 1, 0.15), fontsize=\"large\")\n",
" plt.plot();\n",
"\n",
"# create interactive variables\n",
"N = IntSlider(min=1.0, max=100.0, step=1.0, value=80, readout_format='')\n",
"K = IntSlider(min=1.0, max=N.value, step=1.0, value=40, readout_format='')\n",
"n = IntSlider(min=1.0, max=N.value, step=1.0, value=30, readout_format='')\n",
"\n",
"# enforce K<=N and n<=N\n",
"def update_range(*args):\n",
" K.max = N.value\n",
" n.max = N.value\n",
"N.observe(update_range, 'value')\n",
"\n",
"# display the interactive plot\n",
"interact(hypergeometric_pmf, N=N, K=K, n=n);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Geometric distribution\n",
"\n",
"A random trial has a probability of success equal to $p$ and probability of failure $q=1-p$.Consider the following experiment: we are doing consecutive random trials until we reach a success.The set of outcomes has the form $s=\\overbrace{FF\\cdots F}^{n-1} S$ where number of $F$'s can be any number $n=1,2,\\dots$. Due to independence, $P(s)=q^np$. \n",
"\n",
"Let $X(s)$ denote the number of trials it took to reach success\n",
"\n",
"$$X(\\overbrace{FF\\cdots F}^{n-1} S)=n.$$\n",
"\n",
"Observe that \n",
"\n",
"$$f(n)=q^{n-1}p,\\quad n=1,2,\\dots.$$\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Definition (Geometric distribution)**\n",
"\n",
"
\n",
"\n",
"The pmf of the random variable $X$ is called Geometric distribution.\n",
"\n",
"
\n",
"\n",
"Suppose $X$ is a random variable with mean $\\mu_X$ and standard deviation $\\sigma_Y$. Let $Y=aX+b$ where $a$ and $b$ are any two numbers. Then mean and the standard deviation of $Y$ are\n",
"$$\\mu_Y=a\\mu_X+b,\\quad \\sigma_Y=|a|\\sigma_X.$$\n",
"\n",
"
\n",
"\n",
"In Bernoulli experiment the outcome is either of Type I (success) with probability $p$ or Type II (fail) with probability $q=1-p$. \n",
"\n",
"
\n",
"\n",
"A **Bernoulli trials** are the outcomes of several independent Bernoulli experiments with the same success probabilities.\n",
"\n",
"**Example**\n",
"Consider the coin flip experiment with a biased coin that lands head with probability $p$ and tail with probability $q=1-p$.A single coin flip will be a Bernoulli experiment.$n$ independent coin flips will be a Bernoulli trial."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Binomial distribution\n",
"\n",
"\n",
"**Definition**\n",
"\n",
"
\n",
"\n",
"Let $X$ denote the number of successes in $n$ Bernoulli trials with success probability $p$. Then we say that $X$ has the \\textbf{binomial distribution} with parameters $(n,p)$. \n",
"\n",
"
\n",
" \n",
"If $X$ has binomial distribution with parameters $(n,p)$ then \n",
"* $\\text{range}(X)=\\{0,1,\\dots, n\\}$.\n",
"* $$f(x)={n\\choose x}p^xq^{(n-x)}$$ ${n\\choose x}$ ways there can be $x$ successes: the rest $n-x$ are failures."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "1bd13f24cdd2437381c913d68f14707e",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"interactive(children=(IntSlider(value=50, description='n', min=1, readout_format=''), FloatSlider(value=0.5, d…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# nbi:hide_in\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from ipywidgets import interact, FloatSlider, IntSlider\n",
"from scipy.special import comb\n",
"\n",
"def binomial_pmf(n, p):\n",
" q = 1-p\n",
" range_x = np.arange(0, n, 1)\n",
"\n",
" def binom_pmf(x):\n",
" pmf_val = comb(n, x, exact=True) * p**x * q**(n-x)\n",
" return pmf_val\n",
" mean = n*p\n",
"\n",
" pmf_values = np.array([binom_pmf(x) for x in range_x])\n",
" \n",
" # plot setup\n",
" plt.figure(figsize=(14,7)) \n",
" plt.axhline(y=0, color='k')\n",
" plt.ylim(-0.02, max(np.max(pmf_values)+0.05, 0.2))\n",
" plt.xlim(-1, n+1)\n",
" plt.xticks(np.arange(0, n+1, 5))\n",
"\n",
" # PLotting with plt.bar instead of plt.hist works better when f(x) are knowwn\n",
" plt.scatter(np.array([mean]),np.zeros(1), color =\"blue\", s=200)\n",
" plt.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20)\n",
" plt.bar(range_x, pmf_values, width=1, color=(0.1, 0.1, 1, 0.1), edgecolor='blue', linewidth=1.3)\n",
" #plt.title(\"Binomial distribution\")\n",
" plt.figtext(0.8,0.8, \" n={} \\n p={}\".format(n, p), ha=\"left\", va=\"top\",\n",
" backgroundcolor=(0.1, 0.1, 1, 0.15), fontsize=\"large\")\n",
" plt.plot();\n",
"\n",
"# create interactive variables\n",
"p = FloatSlider(min=0.01, max=1, step=0.01, value=0.5, readout_format='')\n",
"n = IntSlider(min=1, max=100, step=1, value=50, readout_format='')\n",
"# display the interactive plot\n",
"interact(binomial_pmf, n=n, p=p);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Binomial vs hypergeometric\n",
"\n",
"In the hypergeometric distribution with paramaters $(N,K,n)$ when $N>>n$, it is very close to the binomial distribution with paameters $(n,\\frac{K}{N})$. In the plot below you can increase the value of the $N$ while leaving $n$ fixed to see that their histograms approach one another. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "0c11ebd4e4a84e9e9c0512f5f4739f0f",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"interactive(children=(IntSlider(value=200, description='N', max=200, min=1, readout_format=''), IntSlider(valu…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# nbi:hide_in\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from scipy.special import comb\n",
"from ipywidgets import interact, IntSlider\n",
"\n",
"def hypergeo_binom_pmf(N=80,K=40,n=30):\n",
" p = K/N\n",
" q = 1-p\n",
" brange_x = np.arange(0, n, 1)\n",
" hrange_x = np.arange(max(0, n-(N-K)), min(n, K)+1, 1)\n",
"\n",
" def hyper_pmf(N,K,n,i):\n",
" pmf_val = comb(K, i, exact=True) * comb(N-K, n-i, exact=True) / comb(N, n, exact=True)\n",
" return pmf_val\n",
" mean = n * K / N\n",
"\n",
" hpmf_values = np.array([hyper_pmf(N,K,n,i) for i in hrange_x])\n",
" \n",
" def binom_pmf(x):\n",
" pmf_val = comb(n, x, exact=True) * p**x * q**(n-x)\n",
" return pmf_val\n",
" mean = n*p\n",
"\n",
" bpmf_values = np.array([binom_pmf(x) for x in brange_x])\n",
"\n",
" # plot setup\n",
" plt.figure(figsize=(14,7)) \n",
" plt.axhline(y=0, color='k')\n",
" plt.ylim(-0.02, max(np.max(hpmf_values)+0.05, np.max(bpmf_values)+0.05, 0.2))\n",
" plt.xlim(-2,max(np.max(hrange_x), np.max(brange_x))+2)\n",
"\n",
" # Plotting hypergeometric\n",
" plt.bar(hrange_x, hpmf_values, width=1, color=(0.1, 0.1, 1, 0.1), edgecolor='blue', \n",
" linewidth=1.3, label=\"Hypergeometric\")\n",
" \n",
" # Plotting binomial\n",
" plt.bar(brange_x, bpmf_values, width=1, color=(0.3, 0.5, 0.3, 0.1), edgecolor='green',\n",
" linewidth=1.3, label=\"Binomial\")\n",
" plt.figtext(0.78,0.75, \" N={} \\n K={} \\n n={}\".format(N,K,n), ha=\"left\", va=\"top\",\n",
" backgroundcolor=(0.1, 0.1, 1, 0.15), fontsize=\"large\")\n",
" plt.legend()\n",
" plt.plot();\n",
"\n",
"# create interactive variables\n",
"N = IntSlider(min=1.0, max=200.0, step=1.0, value=300, readout_format='')\n",
"K = IntSlider(min=1.0, max=N.value, step=1.0, value=50, readout_format='')\n",
"n = IntSlider(min=1.0, max=N.value, step=1.0, value=20, readout_format='')\n",
"\n",
"# enforce K<=N and n<=N\n",
"def update_range(*args):\n",
" K.max = N.value\n",
" n.max = N.value\n",
"N.observe(update_range, 'value')\n",
"\n",
"# display the interactive plot\n",
"interact(hypergeo_binom_pmf, N=N, K=K, n=n);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### cdf of the binomial distribution\n",
"\n",
"The cdf of the binomial distribution can be computed as \n",
"$$F(x)=E[X\\leq x]=\\sum_{y=0}^{\\lfloor x\\rfloor}{n\\choose x}p^xq^{(n-x)}\\text{ for } x\\in [0,n].$$\n",
"\n",
"However, $F(x)$ does not have a simple closed-form formula. In old times people would use statistical tables that show precomputed values of the cdf for different choices of parameters $(n,p)$. Nowdays, computers can give you very precise values for the cdf for any $n,p,x$.\n",
"\n",
"In the plot below, the shaded area is the value of the cdf at point $x$. The point $x$ is the black dot on the plot. The value of the cdf up to 4 digit precision is given at the top of the plot. "
]
},
{
"cell_type": "code",
"execution_count": 103,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "15c40ff3803a4934a4c163e4088549b3",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"interactive(children=(IntSlider(value=3, description='n', min=1, readout_format=''), FloatSlider(value=0.25, d…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# nbi:hide_in\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from ipywidgets import interact, FloatSlider, IntSlider\n",
"from scipy.special import comb\n",
"from scipy.stats import binom\n",
"\n",
"def binomial_cdf(n, p, x):\n",
" q = 1-p\n",
" range_x = np.arange(0, n, 1)\n",
"\n",
" def binom_pmf(x):\n",
" pmf_val = comb(n, x, exact=True) * p**x * q**(n-x)\n",
" return pmf_val\n",
" mean = n*p\n",
"\n",
" pmf_values = np.array([binom_pmf(x) for x in range_x])\n",
" \n",
" # plot setup\n",
" plt.figure(figsize=(20,10)) \n",
" plt.axhline(y=0, color='k')\n",
" plt.ylim(-0.02, max(np.max(pmf_values)+0.05, 0.2))\n",
" plt.xlim(-4, n+4)\n",
"\n",
" # PLotting with plt.bar instead of plt.hist works better when f(x) are knowwn\n",
" plt.scatter(x,np.zeros(1), color =\"black\", s=100)\n",
" plt.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20)\n",
" barlist = plt.bar(range_x, pmf_values, width=1, color=(1, 1, 1, 0), edgecolor='blue', linewidth=1.3, label=\"Histogran\")\n",
" for ind, val in enumerate(range_x):\n",
" if val<=x:\n",
" barlist[ind].set_color((0.1, 0.1, 1, 0.15))\n",
" barlist[ind].set_edgecolor(\"blue\") \n",
" cdf = binom.cdf(x, n, p)\n",
" plt.title(\"cdf and pdf of binomial distribution\" )\n",
" plt.figtext(0.8,0.8, \" n={} \\n p={} \\n x={} \\n cdf = {:.4f}\".format(n, p, x, cdf), ha=\"left\", va=\"top\",\n",
" backgroundcolor=(0.1, 0.1, 1, 0.15), fontsize=\"large\")\n",
" plt.plot();\n",
"\n",
"# create interactive variables\n",
"p = FloatSlider(min=0.01, max=1, step=0.01, value=0.25, readout_format='')\n",
"n = IntSlider(min=1, max=100, step=1, value=3, readout_format='')\n",
"x = FloatSlider(min=-4, max=10, step=0.1, value=1.75, readout_format='')\n",
"\n",
"# enforce values for x\n",
"def update_range(*args):\n",
" x.max = max(n.value+4,10)\n",
"n.observe(update_range, 'value')\n",
"\n",
"# display the interactive plot\n",
"interact(binomial_cdf, n=n, p=p, x=x);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Negative binomial distribution\n",
"\n",
"\n",
"**Definition**\n",
"\n",
"
\n",
"\n",
"Let $X$ be the number of independent Bernoulli trials with success probability $p$ needed until exactly r successes occur. Then we say $X$ has the **negative binomial distribution** with parameters $(r,p)$.\n",
"
\n",
"\n",
"\n",
"* $\\text{range}(X)=\\{r,r+1,\\dots\\}.$\n",
"*If $X(s)=x$, means there are $r-1$ successes and $x-r$ fails in the first $x-1$ trials and the $x-th $ trial is a success. Therefore\n",
"\n",
"$$f(x)={x-1\\choose r-1}p^{r}q^{x-r}.$$\n",
"\n",
"* When $r=1$, the negative binomial distribution is the same as the geometric distribution. "
]
},
{
"cell_type": "code",
"execution_count": 166,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "4f8724ce4b254192b9fc7fac1afbfdeb",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"interactive(children=(IntSlider(value=30, description='r', max=50, min=1, readout_format=''), FloatSlider(valu…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# nbi:hide_in\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from ipywidgets import interact, FloatSlider\n",
"\n",
"\n",
"def negbinomial_pmf(r, p):\n",
" q = 1-p\n",
" N=100\n",
" range_x = np.arange(r, N, 1)\n",
"\n",
" def negbin_pmf(x):\n",
" pmf_val = comb(x-1, r-1, exact=True) * p**r * q**(x-r)\n",
" return pmf_val\n",
" mean = r/p\n",
"\n",
" pmf_values = np.array([negbin_pmf(x) for x in range_x])\n",
"\n",
" # plot setup\n",
" plt.figure(figsize=(14,7)) \n",
" plt.axhline(y=0, color='k')\n",
" plt.ylim(-0.02, max(np.max(pmf_values)+0.05, 0.1))\n",
" plt.xlim(np.min(range_x)-2, N+1)\n",
"\n",
" # PLotting with plt.bar instead of plt.hist works better when f(x) are knowwn\n",
" plt.scatter(np.array([mean]),np.zeros(1), color =\"blue\", s=200)\n",
" plt.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20)\n",
" plt.bar(range_x, pmf_values, width=1, color=(0.1, 0.1, 1, 0.1), edgecolor='blue', linewidth=1.3)\n",
" plt.title(\"Negative binomial distribution\")\n",
" plt.figtext(0.8,0.8, \" r={} \\n p={} \".format(r,p), ha=\"left\", va=\"top\",\n",
" backgroundcolor=(0.1, 0.1, 1, 0.15), fontsize=\"large\")\n",
" plt.plot();\n",
"\n",
"# create interactive variables\n",
"r = IntSlider(min=1, max=50, step=1, value=30, readout_format='')\n",
"p = FloatSlider(min=0.01, max=1, step=0.01, value=0.5, readout_format='')\n",
"\n",
"# display the interactive plot\n",
"interact(negbinomial_pmf, r=r, p=p);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Poisson distribution\n",
"\n",
"\n",
"**Definition**\n",
"\n",
"
\n",
"\n",
"A discrete random variable has Poisson distribution with parameter $\\lambda>0$ if \n",
"\n",
"$$\\text{range}(X)=\\{0,1,2,3,\\dots\\}$$\n",
"\n",
"and \n",
"\n",
"$$f(x)=\\frac{e^{-\\lambda}\\lambda^x}{x!},\\quad x\\in\\text{range}(X).$$\n",
"\n",
"$\\lambda$ is called the **expected rate**.\n",
"
\n",
"\n",
" Notice that $f(x)$ is indeed a pmf:\n",
" \n",
"$$\\sum_{x\\in \\text{range}(X)}f(x)=\\sum_{x=0}^\\infty \\frac{e^{-\\lambda}\\lambda^x}{x!}=e^{-\\lambda} \\sum_{x=0}^\\infty \\frac{\\lambda^x}{x!}=e^{-\\lambda}e^{\\lambda}=1.$$\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "fbe79733d35b4d08b521e57f9bcff687",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"interactive(children=(FloatSlider(value=30.0, description='$\\\\lambda$', readout_format='', step=0.5), Output()…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# nbi:hide_in\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from scipy.special import factorial\n",
"from ipywidgets import interact, FloatSlider\n",
"\n",
"\n",
"def poisson_pmf(lam):\n",
" N=50\n",
" range_x = np.arange(0, N, 1)\n",
"\n",
" def poiss_pmf(x):\n",
" pmf_val = np.exp(-lam) * np.power(lam, x) / factorial(x)\n",
" return pmf_val\n",
" mean = lam\n",
"\n",
" pmf_values = np.array([poiss_pmf(x) for x in range_x])\n",
"\n",
" # plot setup\n",
" plt.figure(figsize=(14,7)) \n",
" plt.axhline(y=0, color='k')\n",
" plt.ylim(-0.02, max(np.max(pmf_values)+0.05, 0.1))\n",
" plt.xlim(0, N+1)\n",
" plt.xticks(np.arange(0, N+1, 5))\n",
"\n",
" # PLotting with plt.bar instead of plt.hist works better when f(x) are knowwn\n",
" plt.scatter(np.array([mean]),np.zeros(1), color =\"blue\", s=200)\n",
" plt.scatter(range_x,np.zeros(range_x.shape), color =\"red\", s=20)\n",
" plt.bar(range_x, pmf_values, width=1, color=(0.1, 0.1, 1, 0.1), edgecolor='blue', linewidth=1.3)\n",
" plt.title(\"Poisson distribution\")\n",
" plt.figtext(0.8,0.8, r\"$\\lambda$={}\".format(lam), ha=\"left\", va=\"top\",\n",
" backgroundcolor=(0.1, 0.1, 1, 0.15), fontsize=\"large\")\n",
" plt.plot();\n",
"\n",
"# create interactive variables\n",
"lam = FloatSlider(min=0.0, max=100, step=0.5, value=30, readout_format='', description=r\"$\\lambda$\")\n",
"\n",
"# display the interactive plot\n",
"interact(poisson_pmf, lam=lam);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Motivation \n",
"\n",
"The Poisson distribution arises in the following way (made up scenario): grocery route bus can arrive every $\\frac{1}{n}$ hour. It is a very unpredictable bus but we know that there are on average $\\lambda$ buses arriving every hour. Assuming whether a bus arrives or not at given time is independent from other times, this becomes a binomial distribution with parameter $p=\\frac{\\lambda}{n}$. We know that the expect number of successes in a binomial distribution with parameters $(n,p)$ is $np$ and so $p=\\frac{\\lambda}{n}$. Therefore, the probability of $x$ number of buses arriving in an hour will be\n",
"\n",
"$$f(x)={n\\choose x}p^x(1-p)^{n-x}={n\\choose x}\\left(\\frac{\\lambda}{n}\\right)^x\\left(1-\\frac{\\lambda}{n}\\right)^{n-x},\\quad x=0,1,2,\\dots, n. $$\n",
"\n",
"* Now imagine the bus starts arriving at more irregular times. To accommodate that we can increase the $n$. \n",
"* When $n\\to \\infty$ that will correspond to the buses potentially arriving at any moment in the continuous interval. \n",
"* In that case\n",
"\n",
"$$f(x)=\\lim_{n\\to \\infty}{n\\choose x}\\left(\\frac{\\lambda}{n}\\right)^x\\left(1-\\frac{\\lambda}{n}\\right)^{n-x}, \\quad x=0,1,2,\\dots.$$\n",
"\n",
"To compute this limit notice that \n",
"\n",
"$$\\lim_{n\\to \\infty}{n\\choose x}\\left(\\frac{\\lambda}{n}\\right)^x=\\frac{\\lambda^x}{x!}\\lim_{n\\to \\infty}\\frac{n(n-1)\\cdots (n-x+1)}{•n^x}=\\frac{\\lambda^x}{x!}$$\n",
"\n",
"and (you should know this from calculus)\n",
"\n",
"$$\\lim_{n\\to \\infty}\\left(1-\\frac{\\lambda}{n}\\right)^{n-x}=e^{-\\lambda}.$$\n",
"\n",
"Consequently, if the buses can potentially arrive any moment within the hour, with $\\lambda$ buses arriving in average at every hour, then the probability that $x$ buses will arrive in a given hour is \n",
"$$f(x)=\\frac{e^{-\\lambda}\\lambda^x}{x!},\\quad x=0,1,2,\\dots.$$\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "21b928516720454293bcea5371e0eb83",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"interactive(children=(IntSlider(value=100, description='n', max=200, min=6, readout_format=''), FloatSlider(va…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# nbi:hide_in\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from scipy.special import comb, factorial\n",
"from ipywidgets import interact, IntSlider, FloatSlider\n",
"\n",
"def poiss_binom_pmf(n, lam):\n",
" p = lam/n\n",
" q = 1-p\n",
" N = 50\n",
" brange_x = np.arange(0, n, 1)\n",
" prange_x = np.arange(0, N, 1)\n",
"\n",
" def poiss_pmf(x):\n",
" pmf_val = np.exp(-lam) * np.power(lam, x) / factorial(x)\n",
" return pmf_val\n",
" mean = lam\n",
"\n",
" ppmf_values = np.array([poiss_pmf(x) for x in prange_x])\n",
" \n",
" def binom_pmf(x):\n",
" pmf_val = comb(n, x, exact=True) * p**x * q**(n-x)\n",
" return pmf_val\n",
" mean = n*p\n",
"\n",
" bpmf_values = np.array([binom_pmf(x) for x in brange_x])\n",
"\n",
" # plot setup\n",
" plt.figure(figsize=(14,7)) \n",
" plt.axhline(y=0, color='k')\n",
" plt.xlim(-2,20)\n",
"\n",
" # Plotting hypergeometric\n",
" plt.bar(prange_x, ppmf_values, width=1, color=(0.1, 0.1, 1, 0.1), edgecolor='blue', \n",
" linewidth=1.3, label=\"Poisson\")\n",
" \n",
" # Plotting binomial\n",
" plt.bar(brange_x, bpmf_values, width=1, color=(0.3, 0.5, 0.3, 0.1), edgecolor='green',\n",
" linewidth=1.3, label=\"Binomial\")\n",
" plt.figtext(0.83,0.75, r\" $\\lambda$={}\".format(lam)+\"\\n n = {}\".format(n), ha=\"left\", va=\"top\",\n",
" backgroundcolor=(0.1, 0.1, 1, 0.15), fontsize=\"large\")\n",
" plt.legend()\n",
" plt.plot();\n",
"\n",
"# create interactive variables\n",
"lam = FloatSlider(min=0.0, max=100, step=0.5, value=5, readout_format='', description=r\"$\\lambda$\")\n",
"n = IntSlider(min=np.floor(lam.value)+1, max=200, step=1, value=100, readout_format='')\n",
"\n",
"# enforce K<=N and n<=N\n",
"def update_range(*args):\n",
" n.min = np.floor(lam.value)+1\n",
" \n",
"lam.observe(update_range, 'value')\n",
"\n",
"# display the interactive plot\n",
"interact(poiss_binom_pmf, lam=lam, n=n);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us see that the Poisson distribution with parameter $\\lambda$ and the binomial distribution with parameters $(n,p=\\frac{\\lambda}{n})$ are close to each other when $n$ is large by looking at their histograms.\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.7"
}
},
"nbformat": 4,
"nbformat_minor": 4
}