Add files via upload
This commit is contained in:
Peter Norvig authored and 
GitHub committed
1 parent
6ed2674a49
commit
c1781aa0da
1 file changed
+433
-94
+433
-94
@@ -4,27 +4,27 @@
|
|||||||
"cell_type": "markdown",
|
"cell_type": "markdown",
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"source": [
|
"source": [
|
||||||
"<div align=\"right\"><i>Peter Norvig<br>March 2019</i></div>\n",
|
"<div align=\"right\"><i>Peter Norvig<br>March 2019<br>Updated May 2021</i></div>\n",
|
||||||
"\n",
|
"\n",
|
||||||
"# Pairing Socks\n",
|
"# Pairing Socks\n",
|
||||||
"\n",
|
"\n",
|
||||||
"[Bram Cohen](https://en.wikipedia.org/wiki/Bram_Cohen) posed a problem that I'll restate thusly:\n",
|
"[**Bram Cohen**](https://en.wikipedia.org/wiki/Bram_Cohen) posed a problem:\n",
|
||||||
"\n",
|
"\n",
|
||||||
"> *You have N pairs of socks, all different, in the dryer. \n",
|
"> *You have N pairs of socks, all different, in the dryer. \n",
|
||||||
"You pull random socks out one-by-one, placing each sock in one of C possible places on the countertop. \n",
|
"You pull random socks out one-by-one, placing each sock in one of C possible spaces on the countertop. \n",
|
||||||
"If the latest sock pairs with another on the countertop, you put the pair away in a drawer. We can ask:*\n",
|
"After putting a sock down on the countertop you compare it to the other socks on the countertop, and if there is a match, you put the pair away in a drawer, freeing up two spaces. We can ask:*\n",
|
||||||
"1. *What's the probability that you will pair up all 2N socks before you run out of the C places?*\n",
|
"1. *Given N pairs of socks and C spaces, what's the probability P(N, C) that you will pair up all the socks before you run out of spaces?*\n",
|
||||||
"2. *For a given N, what is the minimum C to have a 50% chance?*\n",
|
"2. *Given N pairs of socks, what is the minimum C such that P(N, C) > 1/2?*\n",
|
||||||
"3. *A 100% chance?*\n",
|
|
||||||
"\n",
|
"\n",
|
||||||
"The third question is easy to solve analytically: to be 100% certain of pairing up N pairs, we would need N + 1 countertop spaces, because it is possible that the first N socks are all different, but the next one will always pair with one of the others.\n",
|
"Here's how I think about computing *P(N, C)*:\n",
|
||||||
"How about the other questions?\n",
|
|
||||||
"- Can I come up with a mathematical formula? **No.**\n",
|
|
||||||
"- Could I calculate an exact answer? **Only for small N.** Enumerating all permutations of socks would take *O*((2N)!), but since I don't care *which* sock is next, only whether it is a match/non-match, I could get that down to *O*(2<sup>(2N)</sup>), but that's still bad.\n",
|
|
||||||
"- Could I do a simulation? **Yes.** And since we're mostly interested in events with probability around 50%, a simulation with, say, 1000 trials is likely to give very good approximations. If we were investigating events with probability 0.0001, we'd need many more trials.\n",
|
|
||||||
"\n",
|
"\n",
|
||||||
"# Sock Pairing Simulation\n",
|
"- Can I come up with a closed form mathematical formula? **No,** I'm not smart enough.\n",
|
||||||
"\n"
|
"- Could I do a Monte Carlo simulation? **Yes**, but the resulting probability numbers would not be precise.\n",
|
||||||
|
"- Could I do a simple brute-force computation? **Yes**, but only for very small values of *N*.\n",
|
||||||
|
"- Could I define an efficient algorithm? **Yes,** that's what this notebook is about.\n",
|
||||||
|
"\n",
|
||||||
|
"\n",
|
||||||
|
"First get some imports out of the way:\n"
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
{
|
{
|
||||||
@@ -35,63 +35,252 @@
|
|||||||
"source": [
|
"source": [
|
||||||
"%matplotlib inline\n",
|
"%matplotlib inline\n",
|
||||||
"import matplotlib.pyplot as plt\n",
|
"import matplotlib.pyplot as plt\n",
|
||||||
"import random\n",
|
"from functools import lru_cache\n",
|
||||||
"from statistics import mean\n",
|
"from itertools import permutations"
|
||||||
"from functools import lru_cache"
|
|
||||||
]
|
|
||||||
},
|
|
||||||
{
|
|
||||||
"cell_type": "code",
|
|
||||||
"execution_count": 2,
|
|
||||||
"metadata": {},
|
|
||||||
"outputs": [],
|
|
||||||
"source": [
|
|
||||||
"def pair_socks(C:int, socks:list) -> bool:\n",
|
|
||||||
" \"\"\"Can we pair up this sock ordering using a countertop of size C?\"\"\"\n",
|
|
||||||
" countertop = []\n",
|
|
||||||
" while socks and len(countertop) < C:\n",
|
|
||||||
" s = socks.pop()\n",
|
|
||||||
" countertop.remove(s) if s in countertop else countertop.append(s)\n",
|
|
||||||
" return not socks\n",
|
|
||||||
"\n",
|
|
||||||
"def random_socks(N):\n",
|
|
||||||
" \"\"\"A random ordering of N pairs of socks.\"\"\"\n",
|
|
||||||
" socks = 2 * list(range(N))\n",
|
|
||||||
" random.shuffle(socks)\n",
|
|
||||||
" return socks\n",
|
|
||||||
"\n",
|
|
||||||
"@lru_cache(None)\n",
|
|
||||||
"def P(N, C, trials=1000) -> float: \n",
|
|
||||||
" \"\"\"Estimated probability of pairing up all N pairs on a countertop of size C.\"\"\"\n",
|
|
||||||
" return mean(pair_socks(C, random_socks(N)) for _ in range(trials))\n",
|
|
||||||
"\n",
|
|
||||||
"def minimum_C(N, p) -> int: \n",
|
|
||||||
" \"\"\"Minimum countertop size needed to have probability `p` of pairing up all socks.\"\"\"\n",
|
|
||||||
" return next(C for C in range(N // 2, 3*N) if P(N, C) > p)"
|
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
{
|
{
|
||||||
"cell_type": "markdown",
|
"cell_type": "markdown",
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"source": [
|
"source": [
|
||||||
"First, we can use `minimum_C` to find the countertop size necessary to give a 50% chance of pairing up all the socks, for various values of *N*. We see that the relation is very nearly linear, and the slope is a little more than 1/2:"
|
"# Brute-Force Approach\n",
|
||||||
|
"\n",
|
||||||
|
"A straightforward brute-force way to compute *P(N, C)* would be to consider every possible ordering in which we could pull the socks from the dryer, and count how many of these orderings are pairable. For example, with *N* = 2 pairs of socks, there are 24 possible orderings:"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 2,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"text/plain": [
|
||||||
|
"[('▢', '▢', '▩', '▩'),\n",
|
||||||
|
" ('▢', '▢', '▩', '▩'),\n",
|
||||||
|
" ('▢', '▩', '▢', '▩'),\n",
|
||||||
|
" ('▢', '▩', '▩', '▢'),\n",
|
||||||
|
" ('▢', '▩', '▢', '▩'),\n",
|
||||||
|
" ('▢', '▩', '▩', '▢'),\n",
|
||||||
|
" ('▢', '▢', '▩', '▩'),\n",
|
||||||
|
" ('▢', '▢', '▩', '▩'),\n",
|
||||||
|
" ('▢', '▩', '▢', '▩'),\n",
|
||||||
|
" ('▢', '▩', '▩', '▢'),\n",
|
||||||
|
" ('▢', '▩', '▢', '▩'),\n",
|
||||||
|
" ('▢', '▩', '▩', '▢'),\n",
|
||||||
|
" ('▩', '▢', '▢', '▩'),\n",
|
||||||
|
" ('▩', '▢', '▩', '▢'),\n",
|
||||||
|
" ('▩', '▢', '▢', '▩'),\n",
|
||||||
|
" ('▩', '▢', '▩', '▢'),\n",
|
||||||
|
" ('▩', '▩', '▢', '▢'),\n",
|
||||||
|
" ('▩', '▩', '▢', '▢'),\n",
|
||||||
|
" ('▩', '▢', '▢', '▩'),\n",
|
||||||
|
" ('▩', '▢', '▩', '▢'),\n",
|
||||||
|
" ('▩', '▢', '▢', '▩'),\n",
|
||||||
|
" ('▩', '▢', '▩', '▢'),\n",
|
||||||
|
" ('▩', '▩', '▢', '▢'),\n",
|
||||||
|
" ('▩', '▩', '▢', '▢')]"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"execution_count": 2,
|
||||||
|
"metadata": {},
|
||||||
|
"output_type": "execute_result"
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"socks = ('▢', '▢', '▩', '▩')\n",
|
||||||
|
"list(permutations(socks))"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "markdown",
|
||||||
|
"metadata": {},
|
||||||
|
"source": [
|
||||||
|
"Of these 24 orderings, there are just 6 *distinct* ones:"
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
{
|
{
|
||||||
"cell_type": "code",
|
"cell_type": "code",
|
||||||
"execution_count": 3,
|
"execution_count": 3,
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"text/plain": [
|
||||||
|
"{('▢', '▢', '▩', '▩'),\n",
|
||||||
|
" ('▢', '▩', '▢', '▩'),\n",
|
||||||
|
" ('▢', '▩', '▩', '▢'),\n",
|
||||||
|
" ('▩', '▢', '▢', '▩'),\n",
|
||||||
|
" ('▩', '▢', '▩', '▢'),\n",
|
||||||
|
" ('▩', '▩', '▢', '▢')}"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"execution_count": 3,
|
||||||
|
"metadata": {},
|
||||||
|
"output_type": "execute_result"
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"set(permutations(socks))"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "markdown",
|
||||||
|
"metadata": {},
|
||||||
|
"source": [
|
||||||
|
"When there are *C* = 2 spaces, the first and last of these 6 orderings can be paired up, but the others cannot, so *P*(2, 2) = 2/6 = 1/3.\n",
|
||||||
|
"\n",
|
||||||
|
"In general there are (2*N*)! orderings, and (2*N*)! / 2<sup><i>N</i></sup> distinct orderings. This is faster-than-exponential growth. With 10 pairs of socks (and I personally have at least that many, although they might not all be in the dryer at once) there are **2 quintillion** permutations, and **2 quadrillion** distinct permutations.\n",
|
||||||
|
"\n",
|
||||||
|
"# Answering Question 1: P(N, C)\n",
|
||||||
|
"\n",
|
||||||
|
"Fortunately, I don't need to look at every possible ordering to answer the questions. If the first three socks I pull are, say, white, black, and red (WBR), then for the purposes of this problem that's exactly the same as pulling RBW or BRW or RBG. The identity of each sock doesn't matter; all that matters is how many socks are on the countertop (here, 3 in each case), and how many pairs remain to be matched. \n",
|
||||||
|
"\n",
|
||||||
|
"Suppose that we have *C* spaces on the countertop (this number never changes), that *k* of these spaces currently hold socks (this number can rise and fall as we go), and that *N* pairs of socks remain to be paired (this number will hopefully be reduced over time). There are four possibile cases:\n",
|
||||||
|
"- We've succeeded (*P(N, C)* = 1) if there are no more socks to pair (*N* = 0) or if there are fewer pairs remaining than the number of spaces (*N* < *C*). When there are more spaces than pairs, the worst that can happen is that there is one of each pair on the counter, with one space left; from there we can always pair all the socks.\n",
|
||||||
|
"- We've failed (*P(N, C)* = 0) if the countertop is full (*k* ≥ *C*); there's nowhere for the next sock from the dryer to go. \n",
|
||||||
|
"- If *k* = 0, then the next sock from the dryer becomes the lone sock on the counter.\n",
|
||||||
|
"- Otherwise we might succceed or fail, depending on the ordering of the remaining socks in the dryer, of which there are *r* = 2*N* - *k*. Two things can happen with the next sock pulled from the dryer; consider both cases and add their weighted probabilities:\n",
|
||||||
|
" - *k* out of *r* remaining socks match a countertop sock, leaving *k* - 1 socks on the counter, and *N* - 1 pairs to match.\n",
|
||||||
|
" - (*r* - *k*) out of *r* remaining socks don't match, giving us *k* + 1 socks on the counter.\n",
|
||||||
|
"\n",
|
||||||
|
"\n",
|
||||||
|
"With that, I can code up a solution and try it:"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 4,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"@lru_cache(None)\n",
|
||||||
|
"def P(N, C, k=0) -> float:\n",
|
||||||
|
" \"\"\"The probability that we can pair up N pairs of socks using C counter spaces,\n",
|
||||||
|
" assuming there are currently k socks on the counter.\"\"\"\n",
|
||||||
|
" if N == 0 or N < C:\n",
|
||||||
|
" return 1\n",
|
||||||
|
" elif k >= C:\n",
|
||||||
|
" return 0\n",
|
||||||
|
" elif k == 0:\n",
|
||||||
|
" return P(N, C, 1)\n",
|
||||||
|
" else:\n",
|
||||||
|
" r = 2 * N - k # Number of remaining socks in the dryer\n",
|
||||||
|
" match = k * P(N - 1, C, k - 1)\n",
|
||||||
|
" nomatch = (r - k) * P(N, C, k + 1)\n",
|
||||||
|
" return (match + nomatch) / r"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 5,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"text/plain": [
|
||||||
|
"0.3333333333333333"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"execution_count": 5,
|
||||||
|
"metadata": {},
|
||||||
|
"output_type": "execute_result"
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"P(2, 2) # The probability of pairing up 2 pairs of socks using 2 spaces"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 6,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"text/plain": [
|
||||||
|
"0.5105559043639538"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"execution_count": 6,
|
||||||
|
"metadata": {},
|
||||||
|
"output_type": "execute_result"
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"P(10, 7) # The probability of pairing up 10 pairs of socks using 7 spaces"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 7,
|
||||||
|
"metadata": {},
|
||||||
"outputs": [
|
"outputs": [
|
||||||
{
|
{
|
||||||
"name": "stdout",
|
"name": "stdout",
|
||||||
"output_type": "stream",
|
"output_type": "stream",
|
||||||
"text": [
|
"text": [
|
||||||
"Slope = 0.52\n"
|
"CPU times: user 24 ms, sys: 1.69 ms, total: 25.7 ms\n",
|
||||||
|
"Wall time: 25.4 ms\n"
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
{
|
{
|
||||||
"data": {
|
"data": {
|
||||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7YAAAHjCAYAAAAXJoP1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xl4VOXZx/HvM5OdbOxL2AxIZJFFAogsgmgRFUTrrn2tG1pRW620avt2fVtrUVvrQsQV64IbIossakkA2QSCgiJLwpawLyGTfTLzvH8kWkoVImHmZA6/z3V5JXMyzH3nuecw/Dwz5xhrLSIiIiIiIiKRyuN0AyIiIiIiIiL1oWArIiIiIiIiEU3BVkRERERERCKagq2IiIiIiIhENAVbERERERERiWgKtiIiIiIiIhLRFGxFREREREQkoinYioiIiIiISERTsBUREREREZGIFuV0A/XRrFkz27Fjx5DWCAaDeDyhz//hqhPuWqWlpTRq1Cjkddy6fpqVah0tXHMCd66fG/cp0Kzqy42zcuNzAvRaFSm13LhPubVWJMxq1apV+621zY93v4gOth07dmTlypUhreHz+UhKSgppjXDWCXet7Oxshg0bFvI6bl0/zUq1jhauOYE718+N+xRoVvXlxlm58TkBeq2KlFpu3KfcWisSZmWM2VaX++mtyCIiIiIiIhLRFGxFREREREQkokXkW5GNMaOB0enp6fh8vpDWKi0tDenjh7tOuGsFAoGQzwjcu36alWodLVxzAneunxv3KdCs6suNs3LjcwL0WhUptdy4T7m1lptmFZHB1lo7E5iZmZl5Wzjeex6u99KHq044a3m9Xq1fhNTSrCKjVjjnBO5bv3DW0awip45bZ+XG30mvVZFRy637lBtruWlWeiuyiIiIiIiIRDQFWxEREREREYloCrYiIiIiIiIS0RRsRUREREREJKIp2IqIiIiIiEhEU7AVERERERGRiKZgKyIiIiIiIhFNwVZEREREREQimoKtiIiIiIiIRDQFWxEREREREYloCrYiIiIiIiIS0aKcbkBERERERETCZ3puIRPnbaCwqJy0Zf9iwsgMxvZJc7qtelGwFREREREROUVMzy3kwWlrKfdX09bso6CoBQ9OWwsQ0eFWb0UWERERERE5Rfx17nr6BVbzdszvmR/zS5pQTLk/wMR5G5xurV50xFZERERERMTl/NUBVsx7jWfK/0bvmHwKbVMerr6WUuIA2FlU7nCH9aNgKyIiIiIi4lJllVUsm/Uybdc9zSC7le2mBb/038a0wBD8R8TBNqnxDnZZfwq2IiIiIiIiLnPQV8byGc/RZdOznEchhd62fJn5VzY2H8mM97/CHwh8c9/4aC8TRmY42G39KdiKiIiIiIi4ROGBYla+P4ne215klNlNQfRp5J/zJOnnXk+ax0s3AG/0v8+KnBqvsyI7xRgzGhidnp6Oz+cLaa3S0tKQPn6464S7ViAQCPmMwL3rp1mp1tHCNSdw5/q5cZ8Czaq+3DgrNz4nQK9VkVLLjftUJNTK232IDfOfZ9De17jU7GdH3OlsGfhrmp01BowHX2nZN/cd0TmZEZ37sWLFCvr37wcQ8bkqIoOttXYmMDMzM/O2pKSkkNcLR41w1glnLa/Xq/WLkFqaVWTUCuecwH3rF846mlXk1HHrrNz4O+m1KjJquXWfaqi1VucVsmH2Uww/8Aa9zSEKEnuw//y/0673JWDMMf+sm2YVkcFWRERERETkVGWtZeG6LWyf/xQXFr/NWaaYgpSz8F34PG27jjhuoHUjBVsREREREZEI4A8EmbtqAwf/9RRjyqdzrimhsNlAKkb9iradhzjdnqMUbEVERERERBqw8qoA7y9ZS8Wip7i8ejbJpoxdrc6l+uJfk9a+v9PtNQgKtiIiIiIiIg1QUVkVb+fkEr3iaa4IziPRVLC33UgSRz1E67TeTrfXoCjYioiIiIiINCC7iyt54cNVNPksixv4mFhTzaH0S0gc9RAtWnR1ur0GScFWRERERESkAdi818ebHy6h41eT+YknmyhjKcn4IfEX/JKmzTo73V6DpmArIiIiIiLioNXbD/HuhwvpueVFfuFdhIkyVHS/htgRE0hp3NHp9iKCgq2IiIiIiEiYWWvJ2biPGR8uYMieKfzBuwQbE4O/100E+o0jsU2G0y1GFAVbERERERGRMKkOBJm9dhfzPv6Ii4pe41HvCgLRsQT6jSdm8D1EJbXE5/M53WbEUbAVEREREREJsQp/gLdXFZCzYB5Xl03lGe8q/LGNsAPuJfqcu6BRU6dbjGgKtiIiIiIiIiFyuNzPq8u2sXrxHH5U9TbPez/DH5dC8JwHiT77dohv7HSLrqBgKyIiIiIichJMzy1k4rwN7Cwqp2VyHF1bJcLWRdxm32W890v8jZpiB/2O6H63QFyy0+26ioKtiIiIiIhIPU3PLeTBaWsp9wcAS0bJcu7c+h79PBvxJ7SAIX8muu+PIaaR0626koKtiIiIiIhIPf1p9nrK/dVc4FnFXVHT6eXJp9A25dGocdx/7x8hOs7pFl0tZMHWGPMicAmw11rbo3ZbE+BNoCOwFbjKWnvIGGOAJ4CLgDLgx9ba1aHqTUREREREpL6stSzevJ+sBZvoV5bD3THT6erZzrZgC37pv41pgSFUV0Zxv0JtyIXyiO3LwFPAK0dsewD42Fr7F2PMA7W3fwmMAk6v/W8AMKn2q4iIiIiISIMSCFrmrNvFsws20mnPPP4vZganxRSQF2zNvVU/YUbwHAJ4AUhLjXe421NDyIKttXahMabjUZsvBYbVfj8FyKYm2F4KvGKttcAyY0yqMaa1tXZXqPoTERERERH5Pir8Ad5dXcCLORs46/CHZMXMIC1mN8EW3VjR/jFuWt6G0qD95v7x0V4mjMxwsONTR7g/Y9vyiLC6G2hZ+30asOOI+xXUblOwFRERERERRxVX+HlxyQ7eWrGV8yrm80bsLFpE78O26g1DH8OTcRH9PR7+lPbvsyK3SY1nwsgMxvZJc7r9U4KpOUgaogevOWI764jP2BZZa1OP+Pkha21jY8ws4C/W2sW12z8GfmmtXfktjzkOGAfQsmXLvlOnTg1Z/wDBYBCPxxPSGuGsE+5aJSUlJCYmhryOW9dPs1Kto4VrTuDO9XPjPgWaVX25cVZufE6AXqsipZab9qmiiiDzt1WzdEcJl9l/MT5mFk3tIQ4nZ7Ctw9UcbHIWGHPS62pW/zZ8+PBV1trM490v3Eds93z9FmNjTGtgb+32QqDdEfdrW7vtv1hrJwOTATIzM+2wYcNC2C74fD6SkpJCWiOcdcJdKzs7m1DPCNy7fpqVah0tXHMCd66fG/cp0Kzqy42zcuNzAvRaFSm13LBPbd1fyrML85m7ahNXm/n8JXYuyYFD0GEIDJ1AymlD6RmCQPs1zer7C3ewnQHcCPyl9uv7R2y/yxgzlZqTRh3W52tFRERERCSc1hYcJisnj8XrNnNT1Hw+iZtHQqAYOp5HWeZdJHQd4XSL8h1CebmfN6g5UVQzY0wB8FtqAu1bxphbgG3AVbV3/4CaS/1spuZyPzeFqi8REREREZGvWWtZkneArJw81m7awp2x83g8YR6xgVLoNAqGToC2fQn4fE63KscQyrMiX/sdP/qv/81Rezbk8aHqRURERERE5EiBoGX+F7uZlJPHzoLt3JMwlxcbzSc6UA5dxtQE2tY9nW5T6ijcb0UWERERERFxTGV1gPdWF/LswnzK9u9gQuJcxiZ8iNf6Md1/CEN+Di26Ot2mfE8KtiIiIiIi4nq+Cj+vL9/OC4u3EFNSwEPJ8xiZ8CGeoMX0vAaG3AdNOzndppwgBVsREREREXGtfb5KXl6yhVeWbqNJZQF/bTyXoXH/wlR7MH1ugMH3QuMOTrcp9aRgKyIiIiIirrP9QBmTF+Xx1soCOgS380LT+fQz/8JUxUD/W+GceyAlzek25SRRsBUREREREdf4YudhsnLymf35Tnp4dvBuszn0OJyDqUiAgXfV/JfU0uk25SRTsBURERERkYhmrWVZ/kEm5eSxcOM+zo7dytyWc+hStAgqkmtOCHX2ndCoqdOtSogo2IqIiIiISEQKWsvcdbvJysljzY4iRjTKI6f1B3Q4tBQqUmH4r6D/OIhPdbpVCTEFWxERERERiShV1UGmrylk0oJNbDlQxtiUPLLazKTVwU+hshmc/zvodyvEJjndqoSJgq2IiIiIiDRY03MLmThvAzuLymmVEseA05qwLP8gu4vLuaHJBt5tM5MmB3MhuhWMfBj6/hhiEpxuW8JMwVZERERERBqk6bmFPDhtLeX+AAC7Dlfw/poCbkj9kgltZpJ8cC1Et4OLH4PeN0B0nMMdi1MUbEVEREREpEF6eM76b0KthyCjPCu4K2o6XSu2Q3xHKn4wkbj+P4aoGGcbFccp2IqIiIiISIOyflcxz+bksae4Ei8BxniWMD7qfTp7drI52IZ7q+7kb3f9EX9ZOXEKtYKCrYiIiIiINADWWj7deohJ2ZtZsGEfKTGW/4ldyM12Oh09e1gfbM+dVfcwN9if1qmNwKsoI/+mZ4OIiIiIiDgmGLR8/NVesnLyWLXtEK0S4OXunzF032t4igtYRzq3Vd3HR8GzsHiIj/YyYWSG021LA6NgKyIiIiIiYecPBHl/zU6ezclj094SOqV6eLvXajILX8Xk7YZ2A2D0E2z2ncGX8zdCUTlpqfFMGJnB2D5pTrcvDUxEBltjzGhgdHp6Oj6fL6S1SktLQ/r44a4T7lqBQCDkMwL3rp9mpVpHC9ecwJ3r58Z9CjSr+nLjrNz4nAC9VkVKrePNqawqwLQ1u5myvIDdxZX0bG6YfuYyeha8hmfDAarbDaRq1BME2p0DxjCiFYw4vd9/PMbXj+/G9QtnLTf9/ReRwdZaOxOYmZmZeVtSUugvuhyOGuGsE85aXq9X6xcLine truncated
|
"text/plain": [
|
||||||
|
"0.22085941648149854"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"execution_count": 7,
|
||||||
|
"metadata": {},
|
||||||
|
"output_type": "execute_result"
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"%time P(400, 200) # The probability of pairing up 400 pairs of socks using 200 spaces"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "markdown",
|
||||||
|
"metadata": {},
|
||||||
|
"source": [
|
||||||
|
"Clearly `P` is fast, and these answers look reasonable. But are they actually correct? Here are some tests with the right answers computed by hand:"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 8,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [],
|
||||||
|
"source": [
|
||||||
|
"for N in range(1, 10):\n",
|
||||||
|
" assert P(N, N + 1) == 1 # N + 1 spaces always works\n",
|
||||||
|
" assert P(N, 1) == 0 # 1 space always fails\n",
|
||||||
|
"assert P(0, 0) == 1 # Zero socks can always be paired\n",
|
||||||
|
"assert P(2, 2) == 1/3 # Only works if 2nd sock matches 1st; there's a 1/3 chance of that\n",
|
||||||
|
"assert P(3, 2) == 1/15 # 1/5 for 2nd sock to match 1st; reduces to f(2, 2) \n",
|
||||||
|
"assert P(4, 2) == P(4, 2, 1) == (1 * P(3, 2, 0) + 6 * P(4, 2, 2)) / 7 # r = 7 remaining socks, only 1 matches"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "markdown",
|
||||||
|
"metadata": {},
|
||||||
|
"source": [
|
||||||
|
"Why is `P` so fast (especially given that the brute-force approach was so slow)? You can think of each call to `P(N, C)` as mostly fetching probabilities from the cache, then perhaps doing a small bit of arithmetic and adding a new entry in the cache. So most of the run time is going to be in cache management. A call to `P(N, C, 0)` will fill cache entries for the first argument ranging from 0 to `N`, and the third ranging from 0 to `C`, so the cache can have no more than *N* × *C* entries, and the run time (and space) requirements should be at worst linear in *N* and *C* (or you could call it quadratic in *N*).\n",
|
||||||
|
"\n",
|
||||||
|
"\n",
|
||||||
|
"# Visualizing P(N, C)\n",
|
||||||
|
"\n",
|
||||||
|
"I can plot *P(N, C)* for various values of the number of pairs of socks, *N*, each denoted by a different color S-shaped curve, and of the number of countertop spaces, *C*, each denoted by a point on a curve:"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 9,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7AAAAHgCAYAAACcrIEcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydeXxU5b3/3+fMmsxk30jYE4ggCAVR3EVQq1YLqHhVbKvt1bYXbUVrtbVX9N6qba2WLlisXfxZl1q0gHpbWotCrVvBVBQIhrAGAiSZbLNktnPO74/JZCOEbWaeWZ63L1/MmWeWT855cnI+5/kuimEYSCQSiUQikUgkEolEkuyoogVIJBKJRCKRSCQSiURyLEgDK5FIJBKJRCKRSCSSlEAaWIlEIpFIJBKJRCKRpATSwEokEolEIpFIJBKJJCWQBlYikUgkEolEIpFIJCmBNLASiUQikUgkEolEIkkJzKIFHC/FxcXGmDFjRMvA6/XicDhEywBA13VUNTnuRUgtg5NMWuTcHRypZXCSSUuyzN1k2idSy+Akk5ZkmbeQXPtFahmcZNIi5+7gSC2DE2stH374YYthGCWDjaWcgR0zZgwbN24ULYN169Yxa9Ys0TIAcLvd5OTkiJYBSC1HIpm0yLk7OFLL4CSTlmSZu8m0T6SWwUkmLckybyG59ovUMjjJpEXO3cGRWgYn1loURdlzpLHksOwSiUQikUgkEolEIpEcBWlgJRKJRCKRSCQSiUSSEkgDK5FIJBKJRCKRSCSSlEAaWIlEIpFIJBKJRCKRpATSwEokEolEIpFIJBKJJCWQBlYikUgkEolEIpFIJCmBNLASiUQikUgkEolEIkkJpIGVSCQSiUQikUgkEklKIA2sRCKRSCQSiUQikUhSAmlgJRKJRCKRSCQSiUSSEkgDK5FIJBKJRCKRSCSSlEAaWIlEIpFIJBKJRCKRpATSwEokEolEIpFIJBKJJCWQBlYikUgkEolEIpFIJClB3Aysoii/VRSlSVGUzUcYVxRF+ZmiKPWKonysKMr0eGnJFOoOubn0J+upO+QWLeVwmmph2VmRf/s+Pp6xJKa+rZ55q+dR31bf73F0bOEbC484NtT74jH2SOMjg46lO4Ht29lx5VUEtm/v9zg6dvA/rj/i2FDvi8fYyWpJFVyNHl586ANcjZ7Dtgcbq/+zfsSxwd53vJ9/LGPpSEvDHp65+79oadjT7/Hxjq144J4Tel/fsUwmdMjLwZ98SOiQt9/jWI1JBqepqYlnnnmGpqamnu1ly5bR1NTU7/HRxiTHj8dTx/sfXIbHU9fv8cAxSWzZ5u3iwg9q2ebtEi3lhFEMw4jPByvKBYAHeNYwjMmDjF8B3AFcAcwEfmoYxsyjfe6MGTOMjRs3xlrucbNu3TpmzZolWgYAbrcbky2LS55YT2OHn4q8LN646wKyrWYhWnJyciKGc8UtsOB3kD8Kls2Ejn2QWxF5YWcj5I2AW9+Cpy86+ljeCFj0AVgdx68lDtS31fOtf3yLH1/wYyqcFcxbPY+D3oOUZZcBcMh3iHJHOS987gVu+L8bkn5s5dyVZFuy47KvhiKexyiwfTv7Ft/FiJ88gWX4cHZceRXhAwcwl0V+9vChQ1jKyxmz4o/sWnBd2oxZysupfP011OzYHM+TPUauRg9/e3oLl946iaIKZ8/27C9NYM2vNuNpC5BTYOfa+2aw4gcb8LQFcObbAPC0DxhrDeAsOMLYgPcd8TOG+vyjjOUU2LlhyUz8QV/c5u3xcqzHp6VhD68v/SFX3nkvQM/jvJIyfnf313G7WsgpLAYM3K0ucotKWPjIEzz33cUJG7v58Sex2O0J3S+JoO/1QuiQF9cL2yi6cQJAz2NTgZ1DT3yI1hFAzbWioKB1BjDl2Si5/TM0/+Kjkxoz5dkou+t0vIHUm7uxoKmpiRUrVrBgwQJKS0t7tufNm8cf//hHOjo6yMvL49Zbb+Xpp5+mo6OD3NxcADo7O486tmjRIqxWa0y0JuvcPVE8njo2b/kGkyf9DKDncVbWCN5//7P4Awew2YYBEAgcxG6v4IwZf2LDhvn4Awew2ys4a+YafD4tafZLMh2jI2nZ5u3iq5t389TkMUxwZPVsL504ils372Z/IMRwu4X1Z07AYTLFVcuJoijKh4ZhzBh0LF4GtvuLxwCvH8HAPgWsMwzjxe7tT4FZhmEcGOozpYHtT90hN1///UaGF2bzwc5WAmEdm1nlklPL+MWNiV/Udrvd5NjUXsOaNwIqPgPb34CwH5TuRX9DB7MNsgqhqxXCgaHHzHY45YqIIT4eLXE4wfhCvh7DWu4o59SiU3l7/9sEtABqd1CDjo5NtZFvz6c90E5AC6B0/6ejY1WtPWNBLSh0zKbauGjURTx24WMx31dHI9bHKGpaKx59hH3fvLPH0NknTcLzj39gBAKgds8zXUex2TAVFKC1tQ0+lp+P1t7ePaYASmTMakUtyEdva8cIBkFRIv9Hx/Lz0dvFjCk2G87Zsxnxkydisk9P5Bgdi0k1m1UMw0ALG5gsKnanBb8nhBbSQQEFMAyGHjMr2J1W/J4gWtjoN3b467o/Y8DrjmfMZFEZO6WYc64fnfQXL30N61Amtayyip3//hAtFJlXkZ/dwGSxkJ2bh6+z87jHsnJy6XK7Bx2zO3PxezrRQqHDxsbNOKvHYMdrvySCgSZ179MfMurW04c0qZbhTvyftkL3HATAAMwKqsOC7g0d/1i2Bd3XPWZWyDq1COtVI5J+7saKYzGpZnPkRn84HMZsNpOdnY3P5yMcDqMokR1qGMZRx0455RQWLFgQE93JZI5O5Fq3r2EdyqTm5EzG5VqHrgfoDQjVUVUbFkshoVAruh5AVW0UF1/MmNEPJ81+SZZjtM3bxa0f7+TpKZVDmtS/nF7N5Rvr2B8IYVcVdAMChoFNVbisKI+nJo+JiZ5MMbCvAz8wDOOf3dtrgXsNwxjSnUoD24svGOaSJ9azv91/2FiWxcRDn5/EdWeMTKgmt9tNzpo74NO/RAwrKpG/pjGYZ5ZsuPwxmH7TsWuJ4S9SdNV1WPYwNh7a2GNYje7/Uhm7yc53Z36X+ePnJ/R7T/YYHWmVVbHZwNAxAsGIKTWMbleTGShZWZR9734KrrnmpD/reI9RKKDxwkPvD25SHRa6PEH0cOoeC7NV5cy5o5g2p1K0FKD/8Yma1s9+/U5e+8mjuF0tQ5rUfo+TALPNxpxbvsbkiy456c8SdYGpB7XDTGq4w485335kkwrdd03iq02xqNg+O5zi88bE94uOkVgfo76rrPn5+SxbtuyoJjVWWCwWLr/8cqZPP/mFg2QxR3Ds17pR0zpxwo/YvHlRz8rpkUwqmIhMeP2YdKhqFqNH3Utl5RdO7AeJMSKP0YmY1CKLGVcoTEA//CSTpSo8Mn4EN1QUnbS2TDGw/wc8OsDAftswjA8Hee1twG0AZWVlp//hD3+Im+ZjxePx4HQ6hWpY9pGfj5o0Qkf4/c+1ws9mH3vI7cmS7d3LlI8ewKJ5MenBuHxH0JLHu+c+e0yv1XUdVY1NmndAD/Bw48O0aW0x+bxkxKk6eXTkown9zpM6RoEAxQ/9D2pbG1pBAeHRo7B9shklHKZ7AS6j0XJyaHnsRyf9OcdyjPwdBvveMRhxrkLzZgN3IxjaSX910mKywYT5yVEDMXp8tFCQrX94hqCnE8VkBgwMTYus1EPK3MAxZ2Uz9eb/OunPieX5/2hY3VC2SeXQVJ2CHQqOJgVVV3pubEaiXQyU7sciCVsNds8+NtMQb2JxjLxeL1u2bGHChAls2bKFQCCAzWYjJycHl8uV0JszFouFc88996Q/J5Fz92gc6VrXMPajG8tRla8BxejG94A2wELEmIY4/C7NyZKDSV0ao886OUQdI78B3yIXFwoWDAwUQt2Pc9HpRCXUczes77lm6KuiXHR+pXSetL5Y75eLLrroiAY28UmSvewD+i4PjgAaB3uhYRi/An4FkRVY0SufIH4F9o8bGtjs2nJE85plMfHfn5/ErBkJWoENemHZ7Rihtvj9ebZkY73iUWZNm3VML4/lnaBvrf8WHiM9C7hAZAX23rPuZda4WQn93pM5RvvuXIzH642EcLW3Y25t7RnLdPOqZGUx4jv3MTkG56ijHaPoimvAHWDfehPhoI6hJccFcjzoWYGdJX4FtqVhD68+8Sifv+s7vPfKi2iBSDSOofVZWUoR4wrdK7BfujUh8/ZkiYYJFy4Yj+u5bWjeAKNrrBA2MPTI/O9rVtUkOCspFhXnpcOZlcIrsIOtsvp8Pj7++OOe1wQCAQKBQKzlDonFYuGKK65g2rRpJ/1ZyboC23+V9b/xBw5gtT1FTs5kWlo83TcL+i5exO7cE1mB/RaVlbNi9pknQyKPUd9c1id2HcLt6sDQDYJ9zikhFFz0zWMdeL458vknS1X4n+rRzCpPvhXYoRBpYF8FblcU5Q9Eijh1HC3/VdLLD9dsoys0+BKHzawyZ2IpCxJlXgFWLwJv0+C/Ioopkt+qhwbkudohuwh8rkHyYweMme1QfRlMW5ion6gnZPiyMZfx5t43CemhQV+noqIqKmEjfFgObIG9gLZA26D5sck0NmvkLOaNmxe/nRlDAtu3s+fLX0Hv7IzkpgLoRzFMqhr5Pxw+vhzYFBxTbDacs2aRf/XVJ7+zhyCa5+oosNLVGQQDAr4TC8k74RzYIcZi8RkDx0wWlTGnFTPujJJY7cYTJuT386cfPIjb1cJLS+4jHApGQoSPE0VVURQFvXu1NhY5sCc6VjX9zJiED8cbPajR8rstaB0Bmp76pDtFAYyukwg7UOjOa49hDmzfMbOCfWIh1qknf5EqimAwyPPPP09HRwfPP/88FRUVeL2R6sqxCAc+mRzY6urqmJjXZEXTfGza9BX8gQPU/PuGbrNq4Pc34vfvP8FPVVEUFcMIc/Qc2NmUlMyNzQ+Line truncated
|
||||||
"text/plain": [
|
"text/plain": [
|
||||||
"<Figure size 1152x576 with 1 Axes>"
|
"<Figure size 1152x576 with 1 Axes>"
|
||||||
]
|
]
|
||||||
@@ -110,34 +299,36 @@
|
|||||||
" plt.grid(b=True, which='major')\n",
|
" plt.grid(b=True, which='major')\n",
|
||||||
" plt.grid(b=True, which='minor', alpha=0.2)\n",
|
" plt.grid(b=True, which='minor', alpha=0.2)\n",
|
||||||
" \n",
|
" \n",
|
||||||
"def plotC(Ns=range(10, 201, 10)):\n",
|
"def plotPs(Ns):\n",
|
||||||
" fig('Number of pairs of socks', 'Size of countertop needed for > 50% paired up')\n",
|
" \"\"\"Plot P(N, C) versus C for various values of N and C.\"\"\"\n",
|
||||||
" X = list(Ns)\n",
|
" fig('C = size of countertop', 'P = probability of pairing up all the socks')\n",
|
||||||
" Y = [minimum_C(N, 0.5) for N in Ns]\n",
|
" for N in Ns:\n",
|
||||||
" plt.plot(X, Y, 'o-')\n",
|
" Cs = range(max(2, N // 2 - 20), min(N + 2, N // 2 + 35))\n",
|
||||||
" plt.plot([X[0], X[-1]], [Y[0], Y[-1]])\n",
|
" Ps = [P(N, C) for C in Cs]\n",
|
||||||
" print('Slope =', round(extent(Y) / extent(X), 2))\n",
|
" plt.plot(Cs, Ps, 'd-', label='N = {}'.format(N))\n",
|
||||||
" \n",
|
" plt.legend()\n",
|
||||||
"def extent(X): return max(X) - min(X)\n",
|
" \n",
|
||||||
" \n",
|
"Ns = [10, 25, *range(50, 401, 50)]\n",
|
||||||
"plotC()"
|
"plotPs(Ns)"
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
{
|
{
|
||||||
"cell_type": "markdown",
|
"cell_type": "markdown",
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"source": [
|
"source": [
|
||||||
"Now, for various numbers of pairs of socks ***N***, plot the probability `P(N, C)` of being able to pair up all the socks:"
|
"For example, the rightmost cyan curve corresponds to *N* = 400, with *C* values ranging from 180 to 235. The cyan point on the *C* = 200 vertical line says that *P*(400, 200) is a little more than 0.2.\n",
|
||||||
|
"\n",
|
||||||
|
"Alternatively, we can have the X-axis show not the countertop size *C*, but rather the ratio *C*/*N*:"
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
{
|
{
|
||||||
"cell_type": "code",
|
"cell_type": "code",
|
||||||
"execution_count": 4,
|
"execution_count": 10,
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"outputs": [
|
"outputs": [
|
||||||
{
|
{
|
||||||
"data": {
|
"data": {
|
||||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7MAAAHjCAYAAADxD0ixAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXl4VOX5sO/ZJ/u+kkAI+2IEVBQQZSnuiFbcLYq1qMXK4oLUahF/baWi1GL9WitatG6IKKS4gCwaWQRZRCQQEggkZJKQbZKZZPbz/TEkksxMksk2k5n3vq65JnnP+57nOXNmOc95NpkkSQgEAoFAIBAIBAKBQNCbkPtaAYFAIBAIBAKBQCAQCLxFGLMCgUAgEAgEAoFAIOh1CGNWIBAIBAKBQCAQCAS9DmHMCgQCgUAgEAgEAoGg1yGMWYFAIBAIBAKBQCAQ9DqEMSsQCAQCgUAgEAgEgl6HMGYFAoFAIBAIBAKBQNDrEMasQCAQCAQCgUAgEAh6HcKYFQgEAoFAIBAIBAJBr0PpawW8JT4+XsrIyOjweqPRSFhYWNcp5AaHw4Fc3r33CYSM9hMI5zxQzoU430JGdyDOefDJEOc8+GSIcx58MsQ5Dz4Z55/zffv2VUiSlNDWml5nzGZkZPD99993eP327duZNGlS1ynkhrq6OiIiIoQMP5ERCOc8UM6FON9CRncgznnwyRDnPPhkiHMefDLEOQ8+Geefc5lMdqo9a0SYsUAgEAgEAoFAIBAIeh3CmBUIBAKBQCAQCAQCQa9DGLMCgUAgEAgEAoFAIOh19LqcWXdYrVaKi4sxmUxtzo2KiiI3N7db9fH3JGytVktaWhoqlaqLtRIIBAKBQCAQCASCniEgjNni4mIiIiLIyMhAJpO1OrcnkpftdjsKhcIvZUiSRGVlJcXFxfTv378bNBMIBAKBQCAQCASC7icgwoxNJhNxcXFtGrICkMlkxMXFtcuLLRAIBAKBQCAQCAT+SkAYs4AwZL1AvFYCgUAgEAgEAoGgtxMwxqxAIBAIBAKBQCAQCIKHoDRmPz1whgkvbKX/UxuZ8MJWPj1wptP7lMlkPPbYY03/L1++nCVLlnR6v08//TTp6emEh4c3Gzebzdx+++0MHDiQSy+9lMLCwk7LEggEAoFAIBAIBILeQtAZsxsPl7F43Y+cqWlAAs7UNLB43Y+dNmg1Gg3r1q2joqKiaxQ9x/Tp09mzZ4/L+JtvvklMTAz5+fksWLCARYsWdalcgUAgEAgEAoFAIPBnAqKa8fk8l/0TR0pqPW4/cLoai11qNtZgtfPk2kO8v+e02zXDUyP54/QRrcpVKpXMmTOHFStWsHTpUu8V98Bll13mdnzDhg0899xzAMycOZNHHnkESZJEPqxAIBAIBAKBQCAICgLOmG2Llobsz+OOTu977ty5ZGVlNQs3bsm2bdtYsGCBy3hoaCg7d+5st6ySkhLS09MBpyEdFRVFZWUl8fHx3isuEAgEAoFAIBAIBL2MgDNm2/KgjvvzV+hqzS7jfaJD+PDBcZ2SHRkZyaxZs1i5ciVhYWFu50yePJmDBw92So5AIBAIBAKBQCAQBDtBlzM7b3IGISpFs7EQlYInrh7SJfufP38+b731Fkaj0e32bdu2MWrUKJfH+PHjvZKTmppKUVERADabDb1eT1xcXKf1FwgEAoFAIBAIBILeQLd5ZmUy2ZvADUC5JEkj3WyXAa8A1wH1wH2SJO3vLn0auX5kElptCC9+eYySmgZSo0N44uoh3DS6T5fsPzY2lpkzZ7Jq1Sruv/9+l+1d5ZmdPn06q1evZty4caxdu5YpU6aIfFmBQCAQCAQCgUAQNHRnmPF/gFeBtz1svxYYdO5xKfD/zj13OzeN7tNlxqs7Fi5cyGuvvdYl+3ryySd57733qK+vJy0tjQceeIAlS5Zw//33c9999zFw4EBiY2P54IMPukReMPPpgTPddpND0AEOrYEtS0FfDFFpMPVZyLqtzflX6ovhQDvmd0RGD7DxxEZe2f8KpcZSksOSmTdmHtdnXt+lawJNhs6oI2VtSrtkbF+1FNXra4jW26mJUmCdcxuTfv2sx/n67GzKV/wNm06HMiWFxAXziZo+vVUZ3q4JFBl535Wya30Bhioz4bEaxs0YwOBLk1uV4e2an+c7OLVpRzfLaP9xeEtuzjZyPnibuoqzRMQnMPGOWQybOLl9ayoriIiLb3NNT8joCYwHyqn9spABNXJ0u/cQeXUGYaMT27XGXmNGEa1pc8358w3tmN8RDh06xJYtW9Dr9URFRTF16lSysrJ8KqMndPJXdKXrOVGwHJNZh1aTQuaAx0lJnuFrtbqdj0ur+MsJHWfMVvpoVCzOTOGW5Fhfq9Upus2YlSTpG5lMltHKlBnA25IkScBumUwWLZPJUiRJ0nWXTt2JwWBo+jspKYn6+vou2e9f//pX/vrXv7qMa7VaPvrooy6RIXAasovX/UiD1Q783LIJEAZtV3DOaAz3xjDNfhSsDc7/9UWw4VGw1MOIm1zn//QpfPEU2BqQtWd+izVNMrIfdf7tI4N244mNLNm5BJPdBIDOqGPJziUAHo00b9e4m//HnX+kwdbAVRlXuZWxqXATL+x5AbPd3K413s7vKRm7//MiiSvXorE6/4/V2zH/7X22g1uDVp+dje6ZZ5FMztfKVlKC7hnnPE9GXc2GbEqfdV3jaGgg8pprXObXfvEFZX/+S7vnd2RNV8rwdOx535Wy7d2j2CzOYoqGKjPb3j0K4NEQdLdm63+PYrXYGXiRqyGRv6+cnDXHsVvbN78ja9zNb+s4wDsjMDdnG5tefxWbxfneras4y6bXXwXwcs1KrGYzg8dd7jI/b9e3bFv9OjaLpV3zPa9pXa+O4I2haTxQTs2640hWBzJk2GvM1Kw7DuBxjWF/GfpP8pHOncPGNZLVTugFCS7z6388iz77hMv81mR4y6FDh8jOzsZqdX756PV6srOzAbrMePRWRk/o1FN4a5jqStdz9OjTOBzO33+TuYSjR58GCGiD9uPSKh4/VkSDw1kMt9hs5fFjzpTF3mzQypy2ZDft3GnM/s9DmPH/gBckSfr23P9bgEWSJH3f2j4vvvhi6fvvm0/Jzc1l2LBh7dKprq6OiIiIds3tKHa7HYVC0fZEH8poz2vWE69VT8jYvn07kyZNcrvN7pDI1dVy9xu70TfYXLb3iQ5hx1NT2pTR3cfhd+fCG49mS8MUQBkCVy6C1FFgKDv3KIe6Uuffp3eBw/V89AhR6bDgcKtTuut8XLX2KnRG1/t5UZooHh39KGa7GbPdjMVuwWQ3YbFb+OT4J9TbXG+eqeQq+kf1x2K3NK0z280Yre7z+XsjE36yc9d2ibhaqIyE9ybJ2DHCw/eiJPHaP+zE17luMobIyVy4GMlixmE2I5ktSGYzNWvW4HBzY1KmUqEeNLBpnmQ247Cc+/uc8edvlCZeTEHmjZg1sWjMVQw4sYHk8lZ/bt2uSVOWMGjrFpe5q3+/A0OVa3FFdYiCYRNSsdTbMDfYMNfbMNdbsTTYqK00QfddgnQYmzkXm+lbcNSBPAKl9nKiU0Zx758nuJ3f0tAEUKo1TJvzCIMvnYDNasFmsWC3WrBZrHz0/O8x1lS77EcTGsrIydMwGY2YjUbM9UZMRgNmo5HainLoxuu11oiIT2DOP95qdU57vxPPN04bkankRFzVD21mNI4GG44GG5LJ+Vy75TSS2e66I4UMVWIoks2BZHU4n889sHXN66SI1pDy1Ngu2deKFSvQ6/Uu41FRUW47XED7vaYOh4OamhreeOMNt44UhUJBfHw8Nput2cNsdv28tqXT+fj6+q2RloYpgFyuJSPjESIjRmCxnMVsqXA+m8uxWCrQ6/chSa7XGFpNKhMm5HT1YXTra9VeT2uFxcbkPUc5a3U97jSNiu/Ht15AF3r+nMtksn2SJF3c1ppeYczKZLI5wByApKSki1qG1EZFRTFw4MB26dQThmZP9HvtrIz8/Hy3X6zn43A4kMu7t0ZYd8rYWWLl4zwrlSYHcVo5twxWcVmKkjMGidxKO7lVdo5V2alvw2b6zzXuK1OfT3e/Vv50LhLLvmbIsX+gcPz8Q2iXazg2ZC7lSVcic9jQmCvQmsrRmsoYmP8mSnvbkQp2uQaLOgaLOobI2lzcvbsloGDAr13GBxSs8mp+62tkfD3p01Z17a7z8btTv2v3XCVKVDIVDVKDxzlZIVkoZc55KpkKpUzJ9rrtHuf/MuaXbsfXVa/zao238zuypuzbj3jwMwnteZ9fqwK+HgFJSSPQGkxo60yE1DqftXUmlDY3F8VukORyJJUKmdns8X1lyboASalCUjkfqJRIShWhmzd7XGO4dabLePhHa72a35E1hq8LOTbkLhwKTdOY3G5myLH3CL8yw60M92ss9CvciO3eX2I1gsUonXsGc43b3QAgU4BCfe6hAvm5Z/0pz2uSR7seYekBz9cr7uZ3ZE3x7iPY6jcD5/8wKFGGTuPCe0cgORxY6w2Ya2ux1Okx1+kpO7AXh83qUY43yFUqFGotCo0GpUaLQq1BodFQlXfE45q08a5e0+Kd27ya39aaix5+vBWt2/GdKIGyAdJ3yVFY23/tIiEhc/Nul5AwJoIkl5DkIClwPssh+qTM45qKoa7vh/ijnucXXNP5to3gvED3RGxsLFqttukREhJCbW0tBQUFOBw/y5fL5fTt25eQkBDq6+ubHg0NDc3muSMuLg65XN7scebMGY/z2zIgoXuvSxyOXUh8giRVIpPFIeNm5PLmHUckyQGU4ZBeAAxu99McDRB17nHc4yyFfFXHFfdAd71W30oqXicUy3nvXxUSV2MmGokS5JxBQQly6lqp+StD4n1Z6zYB9My1qMFgIDw8HIDJkyf7vTH7L2C7JEnvn/v/GDCprTBj4ZntGhmB7pltGTYMIJeBVimn/twd4b6xoYzLjGPLine truncated
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7AAAAHgCAYAAACcrIEcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3xTZfv48c/JbNM9oBRaZosgilBQhoPpHqAobhBQHDhAFBRQ3KAi+EPxwa/yyKPgQpbgQhBQRAVZKrJaVgelbehImmae+/dH2rSlLd00wP32FUnPuq9zcprmyr0UIQSSJEmSJEmSJEmS5O80TR2AJEmSJEmSJEmSJNWETGAlSZIkSZIkSZKkM4JMYCVJkiRJkiRJkqQzgkxgJUmSJEmSJEmSpDOCTGAlSZIkSZIkSZKkM4JMYCVJkiRJkiRJkqQzgq6pA6it6Oho0bZt23odQ1VVNBr/zd39OT5/jg2gsLCQoKCgpg6jSv5+/WR8dSfvvfqR8dWdvPfqR8ZXP/58//n7tZPx1Y8/33vg39fPn2MDb3w7duzIEUI0q2z9GZfAtm3blj///LNex7BYLISEhDRQRA3Pn+Pz59gANmzYQP/+/Zs6jCr5+/WT8dWdvPfqR8ZXd/Leqx8ZX/348/3n79dOxlc//nzvgX9fP3+ODbzxhYaGHqlqvf+m3pIkSZIkSZIkSZJUhkxgJUmSJEmSJEmSpDOCTGAlSZIkSZIkSZKkM8IZ1we2Mi6Xi7S0NOx2e422PxM6LjdmfAEBAcTFxaHX6xutDEmSJEmSJEmSpIZ2ViSwaWlphISE0LZtWxRFqXZ7j8eDVqs9DZHVTWPGJ4TAbDaTlpZGu3btGqUMSZIkSZIkSZKkxuC/1ZC1YLfbiYqKqlHyeq5TFIWoqKga11ZLkiRJkiRJkiT5i7MigQVk8loL8lpJkiRJkiRJknQmOmsS2Nraf9zCVXM2sv+4pUGOpygKEydO9P08a9YsXnjhhXofd+rUqcTHxxMcHFxuucPh4PbbbychIYFevXpx+PDhepclSZIkSZIkSZLkz87JBNbmdDPqoy0cyLIy6qOt2Jzueh/TaDSybNkycnJyGiDCUjfeeCNbtmypsHzBggVERESQnJzMhAkTmDx5coOWK0mSJEmSJEmS5G/OyQR28tJ/yLE6EQJyrA4mffVXvY+p0+kYO3Ysc+bMaYAIS/Xu3ZvY2NgKy1euXMnIkSMBuPXWW1m3bh1CiAYtW5IkSZIkSZIkyZ+cFaMQl/Xiqt38m1FQ5fqsAjtHTthQi3M9h1vl27+PMeDN9TQPDah0n/NbhjL9xi7Vlj1u3Di6du3KpEmTqtxm/fr1TJgwocJyk8nE5s2bqy2jRHp6OvHx8YA3eQ4LC8NsNhMdHV3jY0iSJEmSJEmSJJ1JzroEtjqpuUW+5LWEKrzLq0pgayo0NJQRI0Ywd+5cAgMDK91mwIAB7Ny5s17lAJXWtsrBmSRJkiRJkiRJOpuddQlsdTWlX25NZfrXuylyeXzLAvVaXhrShdt6xte7/PHjx5OUlMSoUaMqXd9QNbBxcXGkpqYSFxeH2+0mPz+fyMjIOsctSZIkSZIkSZLk7866BLY6wy+OZ8O+LNbtzcLhVjHqNAzq3LxBkleAyMhIhg8fzoIFCxg9enSF9Q1VA3vTTTfxv//9jz59+vDVV18xcOBAWQMrSZIkSZIkSdJZrdEGcVIU5b+KomQpivJPFesVRVHmKoqSrCjKX4qiJDVWLCd7fdgFRAcbUIDoYCNv3Nq1QY8/ceLEBhuNeNKkScTFxWGz2YiLi/NNzTNmzBjMZjMJCQnMnj2bmTNnNkh5kiRJklRf5gwrn734B+YMa62X1YbjwAFSbrgRx4EDdV6fefsd9dq/MddXJzk3maErh5Kcm9xk+7+W8Vqd9veH2Ou1f8oahv63K8kpa+pc/t0/3t2k5Tf1+denfLL2cPGWxyBrT932r6+sPTCvd9OVXw/7j1sY+v6fDTaVaFNozFGIFwLXnGL9tUBi8WMs8J9GjKUck0HHR6MuITEmmI9GXYzJUP+KaKu19I9vTEwMNputQeaBfeONN0hLS0NVVdLS0nzHDAgIYMmSJSQnJ7Nlyxbat29f77IkSZKks1tO6hEWTnyEnNQjlf58KjVNQF0OD6vf3cWJzEK+efcvXE5PxWWOyped7FQJnmqzcfTBh3CmpJD64EOoNlud1rsPHarX/vVZf/CBMdhTkjn4wP0V1lfH5rLxyLpHOJh3kHHrxmFzNc3+ma7MWu/fIGWveYCDuSmM+3Hs6T93Ww6PbHySgxqVcRufxGarXaVFSfmHLYebtPymPv+6lo+zEBbfhsmWCotv8/58OhWXT/bepim/HkqmEj2YY2uwqUSbQqMlsEKIn4ETp9hkCPCx8PodCFcUpeJ8MY2kY0wIayb0o2NMyOkqUpIkSZLqzXW8kMw523Adr/5DU9kE1WW3s2zmC5jTU1k+80Vs+XnlfnbZ7b79Tk5Ma5OArvt4D0UWFwiwWZz89uXhCst++niPd1m+w7ss385PH5evyaguAcyYMhWP2QxC4DabyZg6rcHWCyHImDKl/Ppnp6A6HKh2O2pRERmTnym3Pn3yM3ishcUPK+mTJuPJyfGuz8khffJkPBaL73H0qYm4crLRCHDmZJE65ZkavPqlnt/8PCfsJxAIzEVmpm+e3qj7CyFQhYpH9eBSXUz7dZpv/5yiHKZtmobdbafIXUSRuwiby1buUegq9D2m/DIFs93s23fKpilYnJYKjwJnQaWPZ3+ehLkoG6FAji2LKT9PqnS7qo7x7C/PVii/qrIq3X/pUMwKCEUhR4Epy26u3f61Kd+RT4E9jwJ7LgVFJyiwmXl26RDMiiguXzBl6RAKrMeLH5kVH5Zj5R7Pbnwac1GZ8jdOosCS4X0UpPseFuux4udp5R7PfnVS+V8NoSA/tcaPZzc85SvfXJTD9F+ehaK8mj+WPYBqzUJBoFqzYNnY2u1f38eyB6AwCxBQmA0rH63V715TeurLXd6pRGm4qUSbgtKYc4cqitIWWC2EuKCSdauBmUKITcU/rwMmCyH+PNUxe/bsKf78s/wme/bsoXPnzjWOy+PxoNVqa7z96XY64qvtNSthsVgICfHfpH/Dhg3079+/qcOokr9fPxlf3cl7r35kfDWjOj0cn70NT74DbZiRmCd78PPmX3z3Xk7qEVa//To3jJ9MWLMYPpr4MBZzDqFRzYhp34GDO7bhcTnR6vWYQsOwFeTjcbnQ6vUk9OzNDeMn43J4+PTF37HmOgiJCODO6b1Y9/EeDv+Vg8elotVraNc1GiGEd5lboNEqRMcHExodyMEd2aieMp8thApK9d+X6wwarri9I50vbQlA2vgJWNevRzgcKAYDgT17EDV6DKqlAMv6DRR89y24ytQeaLUYz+uILjoaZ2oariNHQFVL1ysKSkgIGqMB1VqIKCqqGETJWBJNMK+6y6Ah7JmnUK/vj83tTfqK3EWcsJxA6AQ2t82XGP6d/Tdbj2/FI0prrTWKhrYhbYkMjMQjPL5E0yM8uFW371+36sbqtFLorvgFiE7RoVE0CARCCATepFUg55mXGk+AqjLFnMvN1jOnJrMcvQmufROS7im3uKCgAKMpCKdb9T48qu+546SfK6yvsM5Tft9Kt6n43LutB6dHxe5SK4QeqNfy4k1dGH5xw4wF1FAsFguhoaHbhBA9K1vflAnsN8CMkxLYSUKIbZVsOxZvM2NiYmJ6fP755+XWh4WFkZCQUOO4hBB+PeDR6YgvOTmZ/Pz8Wu+nqioaTWO2PK8fq9VKcHBwU4dRJX+/fjK+upP3Xv3I+E7NYIGYXRrcRkFgroJGVVA1gsLmgpQEC8HBwXhcTv79fCFOawGG4FBMzWLIP3oQ4fHUKDHT6HTEXzaIohNdsGSA8HjzTn0wuKzePLSx6XBwcdEX6A8dQpuVRW3/EgqNBnd8PLq0NBRPxSbJQq+nqHcvAn//A8XlqrBeNRiwDR6Mae1aNE5npesLr7uOoO++ReOoZL3RSOENNxC0ejUah6PCepdBy++D23Hx2hQCnBVfizwTjH3i1N2aNGhQqfzF0KChvbE9GkWDFm2V/26zbcMlKp6/XtHTL6QfSvGV16ABBZTi/0qW/VjwI05R8fwNioFrwq7xbVvC97MC3+Z9W+m+RsXIdeHXle4jVHRum/fhKUTntqF121ikpFFUya+iSVV5NLf6zzXvRoRhq+R32aSqjMsrQKBFKGUeGk25ZQsCi6rcf4S7ue+MhVJ8wt4LWPwc/qfNwq5UfP0ChIZ71Va+31Xhu4YKQil9/plIrrL8OzQdT3Hm3mN8ykHsVPzdCEDLXUr7cttW9nn0c8/eqsvXdjpF+V6fipRKyw/DwEJNv2r3b3t4MTqPvcJytzaAg23vQhXgEd4pMlW1zHMhfM89Jduo4BGi3D7l/lXLrvPuf4f9SwKpWP4JQhnM+7hUgVsFlwruBnzP1Cig04BeAzqNgk4peQ56reJ7rtMUP1fKPC9e/tNRF85KYgo1wNyBQQ0XbANQVZVBgwZVmcA25SjEaUDZdD8OyKhsQyHE/wH/B94a2JNrOfbs2VOrb8ZlDay3D2337t1rvZ+/1EJURdaC1Y+Mr+7kvVc/53J8WVlZLFmyhNtuu43mzZv7lruOF2L+dC+RtyViXrQXj9WB0Vr6YVKjKoTkaFAP7uSCp+/lt6Wf4XF4P1i5bFbyD1tK5wyvwZfVqtvNsT82ow++AFFcgypUcBbU/Fx0Bu8HW3fZT0mqCwUQGr1vkSLcIES5ZRqPg4T9XxKiPYr7ROU9kDShobRdvAjLpk3kzJ2LKCr9IKkEBhL73DTCb7mF3KVLOf7Kq+VqWWuyvlVN11/ck9QXp6N3ln4Qd+oVdt3Tna098mmpiebaFekYy+SIdj0suFJle08zuxUN963xEHDS+q8GBzDz8lcw6UyY9CYLine truncated
|
||||||
"text/plain": [
|
"text/plain": [
|
||||||
"<Figure size 1152x576 with 1 Axes>"
|
"<Figure size 1152x576 with 1 Axes>"
|
||||||
]
|
]
|
||||||
@@ -149,73 +340,221 @@
|
|||||||
}
|
}
|
||||||
],
|
],
|
||||||
"source": [
|
"source": [
|
||||||
"def plotP(Ns=range(10, 101, 10)):\n",
|
"def plotPs2(Ns):\n",
|
||||||
" fig('Size of countertop, C', 'Probability of pairing up all socks, P(N, C)')\n",
|
" \"\"\"Plot P(N, C) versus C/N for various values of N and C.\"\"\"\n",
|
||||||
|
" fig('C / N = ratio of size of countertop to number of pairs of socks', \n",
|
||||||
|
" 'P = probability of pairing up all the socks')\n",
|
||||||
" for N in Ns:\n",
|
" for N in Ns:\n",
|
||||||
" X = range(max(2, N // 2 - 8), min(N + 2, N // 2 + 18))\n",
|
" Cs = range(max(2, N // 2 - 15), min(N + 2, N // 2 + 25))\n",
|
||||||
" plt.plot(X, [P(N, C) for C in X], 'o-', label='N = {}'.format(N))\n",
|
" Xs = [C / N for C in Cs]\n",
|
||||||
|
" Ps = [P(N, C) for C in Cs]\n",
|
||||||
|
" plt.plot(Xs, Ps, 'd-', label='N = {}'.format(N))\n",
|
||||||
" plt.legend()\n",
|
" plt.legend()\n",
|
||||||
" \n",
|
" \n",
|
||||||
"plotP()"
|
"plotPs2(Ns)"
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
{
|
{
|
||||||
"cell_type": "markdown",
|
"cell_type": "markdown",
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"source": [
|
"source": [
|
||||||
"# Exact Probabilities\n",
|
"This shows that as *N* increases, there is a convergence to similar shaped and positioned curves where there is a rapid change from *P(N, C)* ≊ 0 to *P(N, C)* ≊ 1 in the range 0.4 < *C* / *N* < 0.6. (For lower values of *N*, the shape of the curve is not as steep.)"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "markdown",
|
||||||
|
"metadata": {},
|
||||||
|
"source": [
|
||||||
|
"# Answering Question 2: minimum_C(N)\n",
|
||||||
"\n",
|
"\n",
|
||||||
"I said that I could calculate the probability exactly for small values of ***N***. Let's do that. `Pexact` will call `pair_socks` once for each posssible permutation of the socks, and will use `fractions.Fraction` to get an exact rational number as the probability. This gives us a precise answer, and serves as a check on the inexact function `P`."
|
"For question 2, I can try values of *C* until I find one where *P(N, C)* ≥ 50%. If I was worried about efficiency I could do a binary search, but I'm more worried about clarity and correctness, so I'll do a linear search, with a little trick to set a good starting value."
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
{
|
{
|
||||||
"cell_type": "code",
|
"cell_type": "code",
|
||||||
"execution_count": 5,
|
"execution_count": 11,
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"outputs": [],
|
"outputs": [],
|
||||||
"source": [
|
"source": [
|
||||||
"from fractions import Fraction\n",
|
"def minimum_C(N, p=0.5) -> int: \n",
|
||||||
"from itertools import permutations\n",
|
" \"\"\"Minimum countertop size C needed to have P(N, C) ≥ p.\"\"\"\n",
|
||||||
"\n",
|
" start_C = N // 2 if P(N, N // 2) < p else 0\n",
|
||||||
"@lru_cache(None)\n",
|
" return next(C for C in range(start_C, N + 2) if P(N, C) >= p)"
|
||||||
"def Pexact(N, C) -> float: \n",
|
|
||||||
" \"\"\"Exact probability of pairing up all N pairs on a countertop of size C.\"\"\"\n",
|
|
||||||
" return mean(Fraction(pair_socks(C, list(socks))) \n",
|
|
||||||
" for socks in permutations(list(range(N)) * 2))"
|
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
{
|
{
|
||||||
"cell_type": "code",
|
"cell_type": "code",
|
||||||
"execution_count": 6,
|
"execution_count": 12,
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"outputs": [
|
"outputs": [
|
||||||
{
|
{
|
||||||
"data": {
|
"data": {
|
||||||
"text/plain": [
|
"text/plain": [
|
||||||
"{'1000 trials': 0.487,\n",
|
"7"
|
||||||
" '10,0000 trials': 0.464,\n",
|
|
||||||
" 'exact': Fraction(7, 15),\n",
|
|
||||||
" 'exact as decimal': 0.4666666666666667}"
|
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
"execution_count": 6,
|
"execution_count": 12,
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"output_type": "execute_result"
|
"output_type": "execute_result"
|
||||||
}
|
}
|
||||||
],
|
],
|
||||||
"source": [
|
"source": [
|
||||||
"N, C = 5, 4\n",
|
"minimum_C(10) # How many spaces do I need to have at least a 50% chance of pairing up 10 pairs of socks?"
|
||||||
"{'1000 trials': P(N, C, trials=1000), \n",
|
]
|
||||||
" '10,0000 trials': P(N, C, trials=10000), \n",
|
},
|
||||||
" 'exact': Pexact(N, C),\n",
|
{
|
||||||
" 'exact as decimal': float(Pexact(N, C))}"
|
"cell_type": "code",
|
||||||
|
"execution_count": 13,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"text/plain": [
|
||||||
|
"208"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"execution_count": 13,
|
||||||
|
"metadata": {},
|
||||||
|
"output_type": "execute_result"
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"minimum_C(400) # 400 pairs of socks?"
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
{
|
{
|
||||||
"cell_type": "markdown",
|
"cell_type": "markdown",
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"source": [
|
"source": [
|
||||||
"We see that the `P` estimate with 1000 trials is not too bad, and with 10,000 trials it is much better."
|
"# Visualizing minimum_C(N)\n",
|
||||||
|
"\n",
|
||||||
|
"What's the minimum *C* that causes `P(N, C)` to be at least 50% (or in general any probability value)? We can visualize this as a plot:"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 14,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA6wAAAHgCAYAAAChPuB6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdfZhcZ33m+fvnchlKVhEBwj0MOJgKFtUh4c1lIIEEyUCiVtKgJMAkZGdYQ2SHkPDiQNQKBmFmme2eEGCBayYyAS8vs+sQkhg1oRscbDkLw4tlXkw03StnhBlYEzkKBpdaHmyXfvuHquVu6VSdp7ueU+dU1fdzXXWp+zxHOrc7zyXyU51Tt7m7AAAAAAAomnPyDgAAAAAAQBIGVgAAAABAITGwAgAAAAAKiYEVAAAAAFBIDKwAAAAAgEJiYAUAAAAAFNK5eQcIsXnzZr/oootyu/7S0pLOP//83K4f4uTJkzrnnGL/+wMZ42A/xkHGONiPcZAxDvZjHGSMg/0YBxnjGIT9eNtttx1z98eceXwgBtaLLrpIBw8ezO36Bw4c0NatW3O7fohms6lqtZp3jK7IGAf7MQ4yxsF+jIOMcbAf4yBjHOzHOMgYxyDsRzP7TtLxYv9TAAAAAABgZDGwAgAAAAAKiYEVAAAAAFBIDKwAAAAAgEJiYAUAAAAAFBIDKwAAAACgkBhYAQAAAACFxMAKAAAAACgkBlYAAAAAQCExsAIAAAAACuncvAN0Y2aTkiZrtZqazWZuOVqtVq7XD7G0tJR3hFRkjIP9GAcZ42A/xkHGONiPcZAxDvZjHGSMYxD2YyeFHljdfVbSbKPR2FWtVnPLUSqVlOf1Q5ExjqJnZD/GQ8besR/jIWPv2I/xkLF37Md4yNi7QdmPSbglGAAAAABQSAysAAAAAIBCYmAFAAAAgCF1+GhTb/nCCR0+OpjPsDKwAgAAAMAQOnH/g7r8uq/qruOuy6+7VSfufzDvSGvGwAoAAAAAQ+jNn7xdx47fL5d07PiP9UefvD3vSGvGwAoAAAAAQ+YTt35XNy3crR8/eFKS9OMHT+rzC3frE7d+N+dka8PACgAAAABDZmZ+Ufc90Fp17L4HWpqZX8wp0fowsAIAAADAkNm9va5KubTqWKVc0tREPadE68PACgAAAABD5uWXXqjLxi/Qw849NfI97Nxz9ILxC/SyxoU5J1sbBlYAAAAAKKjDR5v6pffc0rGWptv6n7z0qdq88TxJ0uaND9N/fOlTM82aBQZWAAAAACig5VqaO+4+nlhLk7a+4bxzdd3lz9LjNpquu/xSbTjv3H7Gj4KBFQAAAAAK6HQtjSfX0qStS9KWsare+bwN2jJW7VfsqBhYAQAAAKBg0mpphqW2Jk2h3xM2s0lJk7VaTc1m8j3b/dBqtXK9foilpaW8I6QiYxzsxzjIGAf7MQ4yxsF+jIOMcbAf4xjljNNzC4m1NNNzC5qob0pdX2kQ9mMnhR5Y3X1W0myj0dhVreb3FnapVFKe1w9FxjiKnpH9GA8Ze8d+jIeMvWM/xkPG3rEf4xnVjFMT49q7/9CqobRSLmnPjnFVq9XU9ZUGZT8m4ZZgAAAAACiYtFqaYamtScPACgAAAAAFtFxLY0qupUlbHwYMrAAAAACQgV46VKWHamkuHtuYWEuTtj4MGFgBAAAAILJeO1SXbRmr6nNvfH7HWpq09UHHwAoAAAAAkcXoUAUDKwAAAABERYdqPAysAAAAABDRzPxiYkfqzPxi0DoewsAKAAAAABHt3l5XpVxadaxSLmlqoh60jocwsAIAAABARHSoxsPACgAAAAAJQmppdu47mLhOh2ocDKwAAAAAcIbQWpojx04krtOhGgcDKwAAAACcIbiWRp1raUa9QzUGBlYAAAAAWIFamuIo9PvOZjYpabJWq6nZTL5vvB9arVau1w+xtLSUd4RUZIyD/RgHGeNgP8ZBxjjYj3GQMQ72Yxx5ZZyeW0isnZmeW9BEfVPqetEMwn7spNADq7vPSpptNBq7qtX83iYvlUrK8/qhyBhH0TOyH+MhY+/Yj/GQsXfsx3jI2Dv2Yzx5ZJyaGNfe/YdWDaWVckl7doyrWq2mrhfNoOzHJNwSDAAAAAArUEtTHAysAAAAAHAGammKgYEVAAAAwNBJ61BNOye0lqa2eQO1NBliYAUAAAAwVNI6VEPPCamlueHKBrU0GWJgBQAAADBU0jpUQ89B/hhYAQAAAAyNkI5UelQHBwMrAAAAgKExM7+Y2JE6M7+4pnNQDAysAAAAAIbG7u11VcqlVccq5ZKmJuprOgfFwMAKAAAAYGiEdKTSozo4Mh1YzexOM/uWmX3DzA62jz3KzG40szvavz4yywwAAAAABs/ho03t3HewYy1Nt0qakI5UelQHQz/eYd3m7k9390b7+ylJn3f3iyV9vv09AAAAAEh6qHLmyLETiZUzaZU0aR2qoecgf3ncEvwSSR9pf/0RSTtzyAAAAACgoE5Xzii5ciakkiatQzX0HOQr64HVJX3OzG4zsyvax8bc/fuS1P71gowzAAAAABgQaZUzVNKMFnP37P5ws3/t7neZ2QWSbpT0B5L2u/umFefc4+5nPcfaHnCvkKSxsbFLrr/++sxypjl+/Lg2btyY2/VDnDx5UuecU+zP0CJjHOzHOMgYB/sxDjLGwX6Mg4xxsB/X73U3Lene+88+/ojzpPdddn7qer8V9ee40iDsx23btt224jHS0zK9Udvd72r/ereZ/Y2kZ0k6amaPdffvm9ljJd3d4fdeK+laSWo0Gr5169Yso3Z14MAB5Xn9EM1mU9VqsW9lIGMc7Mc4yBgH+zEOMsbBfoyDjHGwH9fv6vO/q737D63qSa2US3rri5+irY0LU9f7rag/x5UGYT92ktk/BZjZ+WZWXf5a0i9J+gdJ+yW9sn3aKyV9KqsMAAAAAAZLWuUMlTSjJcv3rsckfcHMvinpq5L+1t3nJU1LepGZ3SHpRe3vAQAAAAyRbrUzaetplTNU0oyOzAZWdz/i7k9rv57i7u9sH/8Xd3+Bu1/c/vUHWWUAAAAA0H9ptTOhtTS1zRsSK2eopBkdxX46GAAAAMDASaudCa2lueHKRsfKGSppRgMDKwAAAIBoqKVBTAysAAAAAKKZmV9c9Qm+knTfAy3NzC8GrQMrMbACAAAAiGb39roq5dKqY5VySVMT9aB1YCUGVgAAAADRUEuDmBhYAQAAAERFLQ1iYWAFAAAAsEovHapSeu0MtTQIxcAKAAAA4LReO1SXpdXOUEuDEAysAAAAAE6L0aEKxMLACgAAAEASHaoonkLfLG5mk5Ima7Wams3k++P7odVq5Xr9EEtLS3lHSEXGONiPcZAxDvZjHGSMg/0YBxnjGNT9OD23kNiROj23oIn6ptT1fmQsmkHIOAj7sZNCD6zuPitpttFo7KpW87u3vVQqKc/rhyJjHEXPyH6Mh4y9Yz/GQ8besR/jIWPvBnU/Tk2Ma+/+Q6uG0kq5pD07xlWtVlPX+5GxiIqecVD2YxJuCQYAAAAgiQ5VFA8DKwAAADBkQmppdu47mLhOhyqKhIEVAAAAGCKhtTRHjp1IXKdDFUXCwAoAAAAMkeBaGnWupaFDFUXBwAoAAAAMCWppMGwYWAEAAIAhMTO/mFg7MzO/GLQOFA0DKwAAADAkdm+vq1IurTpWKZc0NVEPWgeKhoEVAAAAGBLU0mDYMLACAAAAQ4RaGgwTBlYAAACgQNI6VNPOCa2lqW3eQC0NCo+BFQAAACiItA7V0HNCamluuLJBLQ0Kj4EVAAAAKIi0DtXQc4BhwcAKAAAAFEBIRyo9qhg1hb5h3cwmJU3WajU1m53v4c9aq9XK9fohlpaW8o6QioxxsB/jIGMc7Mc4yBgH+zEOMsaxnv04PbeQ2JE6Pbegifqm4HNCDcLPkYxxDMLfj50UemB191lJs41GY1e1mt/99aVSSXlePxQZ4yh6RvZjPGTsHfsxHjL2jv0YDxl7t579ODUxrr37D60aSCvlkvbsGD/9Z4WcsxZF/zlKZIxhUP5+TMItwQAAAEABhHSk0qOKUcPACgAAAER0+GhTb/nCiY61NN0qaUI6UulRxShhYAUAAAAiWa6cueu4J1bOpFXSpHWohp4DDAsGVgAAACCS05UzSq6cCamkSetQDT0HGAYMrAAAAEAEaZUzVNIAa8fACgAAAEQwM7+YWDkzM78YtA7gbAysAAAAQAS7t9dVKZdWHauUS5qaqAetAzgbAysAAAAQQVrlDJU0wNqtaWA1sxeY2aSZlbMKBAAAAAyq5coZKblyhkoaYG2CB1Yz+1NJL5T0HEmfyiwRAAAAUGDdelSXK2cet9ESK2eopAHWpuPAambvMrOfWHHoJyX9saSr218DAAAAIyWtR1U6VTnzzudt6Fg5QyUNEK7bO6x/I+kvzOwPzKwk6aOSvizpG5Ku7Uc4AAAAoEhCelQBxNNxYHX3L7r7dkk/lDTfPvZsd3+au7+vXwEBAACAIqBHFei/jjfNm9m5kn5Z0lFJvybpKjPbJelqd+/LPyWZ2aSkyVqtpmbz7GcE+qXVauV6/RBLS0t5R0hFxjjYj3GQMQ72YxxkjIP9GAcZO5ueW0jsUZ2eW9BEfdOq4+zHOMgYxyDsx066PeV9g07d/rtB0m+7+yvN7F9LeoeZubvvyjqcu89Kmm00Gruq1fzu8S+VSsrz+qHIGEfRM7If4yFj79iP8ZCxd+zHeMiYbGpiXHv3H1o1tFbKJe3ZMX5WHvZjPGTs3aDsxyTdBtYnuPuvmtl5OvXsqtz9Lkm/Y2ZP70s6AAAAoCBefumFuuWOf9bf/bej+vGDJ+lRBfqg24cuXWtm35D0FUnvXrng7t/INBUAAACQgW6VNCHn0KMK9Fe3D116v7s/3d2f4e4f72coAAAAILaQSpq0c+hRBfqrWw/ruWZ2pZnNmdntZvbN9te/a2blfoYLine truncated
|
||||||
|
"text/plain": [
|
||||||
|
"<Figure size 1152x576 with 1 Axes>"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"metadata": {
|
||||||
|
"needs_background": "light"
|
||||||
|
},
|
||||||
|
"output_type": "display_data"
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"def plotCs(Ns, p=0.5):\n",
|
||||||
|
" \"\"\"Plot minimum_C(N) for various values of N.\"\"\"\n",
|
||||||
|
" fig('N = number of pairs of socks', \n",
|
||||||
|
" f'C = size of countertop needed for P(N, C) ≥ {p:.0%}')\n",
|
||||||
|
" Cs = [minimum_C(N, p) for N in Ns]\n",
|
||||||
|
" plt.plot(Ns, Cs, 'd')\n",
|
||||||
|
" \n",
|
||||||
|
"plotCs(range(1, 101))"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "markdown",
|
||||||
|
"metadata": {},
|
||||||
|
"source": [
|
||||||
|
"The points are almost on a straight line. (I don't fully understand the variation.)\n",
|
||||||
|
"\n",
|
||||||
|
"As we discovered in the previous plot, *C* needs to be a little over half of *N* to get *P(N, C)* > 50%.\n",
|
||||||
|
"\n",
|
||||||
|
"One weakness of my implementation of `P(N, C)` is that it can quickly exceed Python's built-in recursion limit. To get around that, we could (a) only use values of `N` and `C` that keep us below the default recursion limit of 1000; or (b) reimplement `P` to not be recursive; or (c) increase the recursion limit:"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 15,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"name": "stdout",
|
||||||
|
"output_type": "stream",
|
||||||
|
"text": [
|
||||||
|
"CPU times: user 25.3 s, sys: 1.15 s, total: 26.5 s\n",
|
||||||
|
"Wall time: 26.5 s\n"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"data": {
|
||||||
|
"text/plain": [
|
||||||
|
"0.5070227508048063"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
"execution_count": 15,
|
||||||
|
"metadata": {},
|
||||||
|
"output_type": "execute_result"
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"import sys\n",
|
||||||
|
"sys.setrecursionlimit(50_000)\n",
|
||||||
|
"\n",
|
||||||
|
"%time P(10_000, 5_022)"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "markdown",
|
||||||
|
"metadata": {},
|
||||||
|
"source": [
|
||||||
|
"Here is a table of `minimum_C` values reaching up to *N* = 2000:"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "code",
|
||||||
|
"execution_count": 16,
|
||||||
|
"metadata": {},
|
||||||
|
"outputs": [
|
||||||
|
{
|
||||||
|
"name": "stdout",
|
||||||
|
"output_type": "stream",
|
||||||
|
"text": [
|
||||||
|
"N = 10; C = 7; C/N = 0.700; C-N/2 = 2\n",
|
||||||
|
"N = 50; C = 29; C/N = 0.580; C-N/2 = 4\n",
|
||||||
|
"N = 100; C = 55; C/N = 0.550; C-N/2 = 5\n",
|
||||||
|
"N = 200; C = 106; C/N = 0.530; C-N/2 = 6\n",
|
||||||
|
"N = 400; C = 208; C/N = 0.520; C-N/2 = 8\n",
|
||||||
|
"N = 600; C = 309; C/N = 0.515; C-N/2 = 9\n",
|
||||||
|
"N = 800; C = 409; C/N = 0.511; C-N/2 = 9\n",
|
||||||
|
"N = 1000; C = 510; C/N = 0.510; C-N/2 = 10\n",
|
||||||
|
"N = 1200; C = 611; C/N = 0.509; C-N/2 = 11\n",
|
||||||
|
"N = 1400; C = 711; C/N = 0.508; C-N/2 = 11\n",
|
||||||
|
"N = 1600; C = 812; C/N = 0.507; C-N/2 = 12\n",
|
||||||
|
"N = 1800; C = 912; C/N = 0.507; C-N/2 = 12\n",
|
||||||
|
"N = 2000; C = 1013; C/N = 0.506; C-N/2 = 13\n",
|
||||||
|
"CPU times: user 43.3 s, sys: 1.97 s, total: 45.3 s\n",
|
||||||
|
"Wall time: 45.3 s\n"
|
||||||
|
]
|
||||||
|
}
|
||||||
|
],
|
||||||
|
"source": [
|
||||||
|
"%%time\n",
|
||||||
|
"for N in [10, 50, 100, *range(200, 2001, 200)]:\n",
|
||||||
|
" C = minimum_C(N)\n",
|
||||||
|
" print(f'N = {N:4}; C = {C:4}; C/N = {C / N:.3f}; C-N/2 = {C - N // 2:2}')"
|
||||||
|
]
|
||||||
|
},
|
||||||
|
{
|
||||||
|
"cell_type": "markdown",
|
||||||
|
"metadata": {},
|
||||||
|
"source": [
|
||||||
|
"I'm going to stop here, but there are several things you could work on if you're interested:\n",
|
||||||
|
"- Can you develop a closed form mathematical formula for *P(N, C)*?. \n",
|
||||||
|
"- Can you run a Monte Carlo simulation? \n",
|
||||||
|
"- Can you display an animation of the sock-pairing process?\n",
|
||||||
|
"- Can you implement a brute-force enumeration, and use it to verify *P(N, C)* on small values?\n",
|
||||||
|
"- Can you reimplement `P` to not be recursive?\n",
|
||||||
|
"- Can you compute a full probability distribution over the number of pairs of socks paired up?\n",
|
||||||
|
"- What else would you like to explore?"
|
||||||
]
|
]
|
||||||
}
|
}
|
||||||
],
|
],
|
||||||
@@ -235,7 +574,7 @@
|
|||||||
"name": "python",
|
"name": "python",
|
||||||
"nbconvert_exporter": "python",
|
"nbconvert_exporter": "python",
|
||||||
"pygments_lexer": "ipython3",
|
"pygments_lexer": "ipython3",
|
||||||
"version": "3.7.2"
|
"version": "3.7.6"
|
||||||
}
|
}
|
||||||
},
|
},
|
||||||
"nbformat": 4,
|
"nbformat": 4,
|
||||||
|
|||||||
Reference in new issue
Block a user