{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# CS 237 Spring 2019\n",
"# Author: Alina Ene (aene@bu.edu)\n",
"# Used in L20"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# CS 237: Double Sixes \n",
"\n",
"In this notebook, we consider the following experiment:\n",
"\n",
"> Experiment: We repeatedly roll a fair 6-sided die until we get a pair of consecutive sixes for the first time. The rolls are independent.\n",
"\n",
"We are interested in the expectation of the following random variable:\n",
"\n",
"> X = number of rolls we perform"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## To be or not to be (geometric)\n",
"\n",
"What is the distribution of X? On a first glance, it has the flavor of a geometric random variable. One reasoning could be the following:\n",
"\n",
"- Define a trial to be two consecutive rolls\n",
"- A trial is a success if the two consecutive rolls are 66\n",
"- Pr(success) = 1/36, since the die is 6-sided and fair, and the rolls are independent\n",
"- The total number of rolls we perform until 66 is equal to the total number of trials until the first success\n",
"- Therefore X ~ Geometric(1/36) and thus Ex(X) = 36\n",
"\n",
"Think: Do you agree with the above argument?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulation\n",
"\n",
"Let us do the simulation and see what we get. The code below estimates Ex(X)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from numpy.random import randint\n",
"import matplotlib.pyplot as plt\n",
"plt.style.use('seaborn')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# simulation for the double sixes experiment\n",
"# a single experiment\n",
"def single_trial():\n",
" prev_roll = randint(1, 7) # fair 6-sided die roll\n",
" num_rolls = 1\n",
" while True:\n",
" curr_roll = randint(1, 7) \n",
" num_rolls = num_rolls + 1\n",
" if prev_roll == 6 and curr_roll == 6:\n",
" return num_rolls\n",
" prev_roll = curr_roll\n",
" return num_rolls\n",
"\n",
"# perform N trials\n",
"N = 100000\n",
"rolls = []\n",
"trial = [i + 1 for i in range(N)]\n",
"for i in range(N):\n",
" num_rolls = single_trial()\n",
" rolls.append(num_rolls)\n",
"\n",
"avg_rolls = sum(rolls)\n",
"avg_rolls = avg_rolls / N\n",
"plt.bar(trial,rolls)\n",
"plt.xlabel(\"Trial\")\n",
"plt.ylabel(\"Number of rolls\")\n",
"plt.title(\"Ex(# rolls until 66) = \" + str(avg_rolls),fontsize=20)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Where is the mistake?\n",
"\n",
"The simulation says the expectation is about 42, whereas the above argument says it is 36. So there are two possibilities:\n",
"\n",
"- There is a mistake in the code\n",
"- There is a mistake in the above argument\n",
"\n",
"The code seems simple enough and correct enough. So maybe the above argument is worth taking a closer look. Indeed, as we have seen in class, a Geometric(p) random variable counts the number of independent trials until the first success, where p is the probability that a single trial is a success. The above reasoning argues (correctly) that X is the number of trials until the first success, where Pr(success) = 1/36. But it does not argue that the trials are independent!\n",
"\n",
"Are the trials independent? Consider the first two trials: the first trial is rolls 1 and 2, and the second trial is rolls 2 and 3. Since the two trials share the second roll of the die, they are dependent.\n",
"\n",
"Therefore the conclusion of the above argument does not follow."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The math (optional)\n",
"\n",
"The random variable X is not geometric. The trials as defined earlier (trial = two consecutive rolls) are not independent, since two consecutive trials share a roll, but trials that are farther apart are independent. So perhaps we can still use the ideas we developed for the geometric distribution, but apply them with a bit of care. As we now illustrate, this is indeed the case: we will use the ideas from lecture to break down the expectation. The event of interest now is whether the first two rolls were sixes.\n",
"\n",
"Let\n",
"- $S_1$ = first roll is a 6\n",
"- $S_2$ = second roll is a 6\n",
"\n",
"Also, for notational convenience, let $p = \\frac{1}{6}$ and $q = \\frac{5}{6}$.\n",
"\n",
"Since the events $\\overline{S}_2$, $S_2 \\cap S_1$, and $S_2 \\cap \\overline{S}_1$ form a partition of the sample space, the law of total Ex gives us that:\n",
"\n",
"$$\n",
"\\begin{align*}\n",
"\\mathrm{Ex}(X) &= \\mathrm{Ex}(X | \\overline{S}_1) \\Pr(\\overline{S}_1) + \\mathrm{Ex}(X | S_1 \\cap S_2) \\Pr(S_1 \\cap S_2) + \\mathrm{Ex}(X | S_1 \\cap \\overline{S}_2) \\Pr(S_1 \\cap \\overline{S}_2)\\\\\n",
"&= \\mathrm{Ex}(X | \\overline{S}_1) q + \\mathrm{Ex}(X | S_1 \\cap S_2) p^2 + \\mathrm{Ex}(X | S_1 \\cap \\overline{S}_2) pq\n",
"\\end{align*}\n",
"$$\n",
"\n",
"Now we consider each conditional expectation in turn. As in class, the key will be the fact that the process is memoryless about what happened in the past, provided that the past is not so recent. More precisely, we only need to remember whether the previous roll was a 6 or not, since the outcomes of the rolls that were further in the past do not affect the number of rolls we will need from this point onward. With this observation in mind, we have\n",
"\n",
"$$\n",
"\\begin{align*}\n",
"\\mathrm{Ex}(X | \\overline{S}_1) &= 1 + \\mathrm{Ex}(X)\\\\\n",
"\\mathrm{Ex}(X | S_1 \\cap S_2) & = 2\\\\\n",
"\\mathrm{Ex}(X | S_1 \\cap \\overline{S}_2) &= 2 + \\mathrm{Ex}(X)\n",
"\\end{align*}\n",
"$$\n",
"\n",
"The formal argument for the above identities is similar to what we saw in lecture, and it is left as an exercise.\n",
"\n",
"Putting everything together, we obtain:\n",
"$$\n",
"\\mathrm{Ex}(X) = (1 + \\mathrm{Ex}(X)) q + 2 p^2 + (2 + \\mathrm{Ex}(X)) pq\n",
"$$\n",
"By solving for $\\mathrm{Ex}(X)$ in the above equation and using that $q = 1-p$, we obtain:\n",
"$$\n",
"\\mathrm{Ex}(X) = \\frac{1+p}{p^2} = 42$$\n",
"\n",
"Recall that the simulation gave us an empirical expectation of approximately 42."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}