Add files via upload
This commit is contained in:
Peter Norvig authored and 
GitHub committed
1 parent
c3bca4d0ae
commit
3fc0e1049c
1 file changed
+243
@@ -0,0 +1,243 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"<div align=\"right\"><i>Peter Norvig<br>March 2019</i></div>\n",
|
||||
"\n",
|
||||
"# Pairing Socks\n",
|
||||
"\n",
|
||||
"[Bram Cohen](https://en.wikipedia.org/wiki/Bram_Cohen) posed a problem that I'll restate thusly:\n",
|
||||
"\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",
|
||||
"If the latest sock pairs with another on the countertop, you put the pair away in a drawer. 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",
|
||||
"2. *For a given N, what is the minimum C to have a 50% chance?*\n",
|
||||
"3. *A 100% chance?*\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",
|
||||
"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",
|
||||
"# Sock Pairing Simulation\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%matplotlib inline\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"import random\n",
|
||||
"from statistics import mean\n",
|
||||
"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",
|
||||
"metadata": {},
|
||||
"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:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Slope = 0.52\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"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": [
|
||||
"<Figure size 1152x576 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"def fig(xlabel, ylabel):\n",
|
||||
" \"\"\"Set up a figure.\"\"\"\n",
|
||||
" plt.figure(figsize=(16,8)); plt.xlabel(xlabel); plt.ylabel(ylabel)\n",
|
||||
" plt.minorticks_on()\n",
|
||||
" plt.grid(b=True, which='major')\n",
|
||||
" plt.grid(b=True, which='minor', alpha=0.2)\n",
|
||||
" \n",
|
||||
"def plotC(Ns=range(10, 201, 10)):\n",
|
||||
" fig('Number of pairs of socks', 'Size of countertop needed for > 50% paired up')\n",
|
||||
" X = list(Ns)\n",
|
||||
" Y = [minimum_C(N, 0.5) for N in Ns]\n",
|
||||
" plt.plot(X, Y, 'o-')\n",
|
||||
" plt.plot([X[0], X[-1]], [Y[0], Y[-1]])\n",
|
||||
" print('Slope =', round(extent(Y) / extent(X), 2))\n",
|
||||
" \n",
|
||||
"def extent(X): return max(X) - min(X)\n",
|
||||
" \n",
|
||||
"plotC()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"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:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"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
|
||||
"text/plain": [
|
||||
"<Figure size 1152x576 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"def plotP(Ns=range(10, 101, 10)):\n",
|
||||
" fig('Size of countertop, C', 'Probability of pairing up all socks, P(N, C)')\n",
|
||||
" for N in Ns:\n",
|
||||
" X = range(max(2, N // 2 - 8), min(N + 2, N // 2 + 18))\n",
|
||||
" plt.plot(X, [P(N, C) for C in X], 'o-', label='N = {}'.format(N))\n",
|
||||
" plt.legend()\n",
|
||||
" \n",
|
||||
"plotP()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Exact Probabilities\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`."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from fractions import Fraction\n",
|
||||
"from itertools import permutations\n",
|
||||
"\n",
|
||||
"@lru_cache(None)\n",
|
||||
"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",
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"{'1000 trials': 0.487,\n",
|
||||
" '10,0000 trials': 0.464,\n",
|
||||
" 'exact': Fraction(7, 15),\n",
|
||||
" 'exact as decimal': 0.4666666666666667}"
|
||||
]
|
||||
},
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"N, C = 5, 4\n",
|
||||
"{'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": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We see that the `P` estimate with 1000 trials is not too bad, and with 10,000 trials it is much better."
|
||||
]
|
||||
}
|
||||
],
|
||||
"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.2"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 2
|
||||
}
|
||||
Reference in new issue
Block a user