Skip to content

Instantly share code, notes, and snippets.

@qpliu

qpliu/20240809.ipynb

Last active August 12, 2024 14:29
Show Gist options
  • Save qpliu/23e5fabee192f78c629145c2be2251c1 to your computer and use it in GitHub Desktop.
Save qpliu/23e5fabee192f78c629145c2be2251c1 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
20240809.ipynb
{
"cells": [
{
"cell_type": "markdown",
"id": "25e7f5a3-1951-4701-aa55-635b7d0f2e35",
"metadata": {},
"source": [
"[2024-08-09 Fiddler](https://thefiddler.substack.com/p/can-you-hack-gymnastics)\n",
"====================\n"
]
},
{
"cell_type": "markdown",
"id": "a0f0f2b1-0470-4730-8cc4-64eccf136161",
"metadata": {},
"source": [
"[Simulations](20240809.go) comparing 160 million pairs of random scores\n",
"give about 96.17% and, for the extra credit, about 91.67%, for the probability\n",
"that the relative rankings using sums are the same as the relative rankings\n",
"using products.\n",
"\n",
" $ go run 20240809.go\n",
" 153866128/160000000 0.961663\n",
" 146666306/160000000 0.916664\n",
"\n",
"Setting things up"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "b54dfbc1-56f2-4d69-ac42-ba675cf78fcd",
"metadata": {},
"outputs": [],
"source": [
"%display latex\n",
"D_A, X_A, D_B, X_B = var('D_A X_A D_B X_B')\n",
"A_s = D_A + X_A\n",
"A_p = D_A * X_A\n",
"B_s = D_B + X_B\n",
"B_p = D_B * X_B"
]
},
{
"cell_type": "markdown",
"id": "0141aa4f-7c33-4cd3-9405-cb261fac9989",
"metadata": {},
"source": [
"The probability that the relative rankings are the same is\n",
"\n",
"$$ \\int_0^{10} \\frac{dX_A}{10} \\int_0^{10} \\frac{dX_B}{10} (H(A_s-B_s)H(A_p-B_p) + H(B_s-A_s)H(B_p-A_p)) $$\n",
"\n",
"but I can't figure out how to express that without getting `TypeError`s. The integrand is"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "67136757-1131-40a6-8e62-6447b91adcef",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<html>\\(\\displaystyle H\\left(D_{A} X_{A} - D_{B} X_{B}\\right) H\\left(D_{A} - D_{B} + X_{A} - X_{B}\\right) + H\\left(-D_{A} X_{A} + D_{B} X_{B}\\right) H\\left(-D_{A} + D_{B} - X_{A} + X_{B}\\right)\\)</html>"
],
"text/latex": [
"$\\displaystyle H\\left(D_{A} X_{A} - D_{B} X_{B}\\right) H\\left(D_{A} - D_{B} + X_{A} - X_{B}\\right) + H\\left(-D_{A} X_{A} + D_{B} X_{B}\\right) H\\left(-D_{A} + D_{B} - X_{A} + X_{B}\\right)$"
],
"text/plain": [
"heaviside(D_A*X_A - D_B*X_B)*heaviside(D_A - D_B + X_A - X_B) + heaviside(-D_A*X_A + D_B*X_B)*heaviside(-D_A + D_B - X_A + X_B)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"heaviside(A_s-B_s)*heaviside(A_p-B_p) + heaviside(B_s-A_s)*heaviside(B_p-A_p)"
]
},
{
"cell_type": "markdown",
"id": "19f6ec54-4a85-4536-bc96-b9ccc8911fa8",
"metadata": {},
"source": [
"A is ahead of B when adding when $X_B < D_A - D_B + X_A$.\n",
"\n",
"A is ahead of B when multiplying when $X_B < D_AX_A/D_B$.\n",
"\n",
"So the relative rankings are the same when $X_B < \\min(D_A-D_B+X_A,D_AX_A/D_B)$ or $X_B > \\max(D_A-D_B+X_A,D_AX_A/D_B)$."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "b33e27d7-b6e1-4f14-876e-99be5e0a0668",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<html>\\(\\displaystyle \\left[X_{A} = D_{B}\\right]\\)</html>"
],
"text/latex": [
"$\\displaystyle \\left[X_{A} = D_{B}\\right]$"
],
"text/plain": [
"[X_A == D_B]"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"solve(D_A-D_B+X_A == D_A*X_A/D_B, X_A)"
]
},
{
"cell_type": "markdown",
"id": "b2b8ef5d-cf75-4402-8972-4169d9ec6d11",
"metadata": {},
"source": [
"Assuming $D_A > D_B$, then when $X_A < D_B$, the relative rankings are the same when\n",
"$X_B < D_AX_A/D_B$ or $X_B > D_A-D_B+X_A$, and when $X_A > D_B$, the relative rankings\n",
"are them same when $X_B < D_A-D_B+X_A$ or $X_B > D_AX_A/D_B$.\n",
"\n",
"If $X_A < D_B$, the probability that $X_B < D_AX_A/D_B$ is $D_AX_A/D_B\\cdot1/10$, and\n",
"the probability that $X_B > D_A-D_B+X_A$ is 0 when $D_A-D_B+X_A > 10$, or $X_A > 10 - (D_A-D_B)$,\n",
"and $(10 - (D_A-D_B+X_A))/10$ when $X_A < 10-(D_A-D_B)$.\n",
"\n",
"If $X_A > D_B$, the probability that $X_B < D_A-D_B+X_A$ is $\\min(1,(D_A-D_B+X_A)/10)$, and\n",
"the probability that $X_B > D_AX_A/D_B$ is 0 when $D_AX_A/D_B > 10$, or $X_A > 10D_B/D_A$,\n",
"is $\\max(0,(10 - D_AX_A/D_B)/10)$.\n",
"\n",
"The probability that the relative rankings are the same is"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "0b137e24-0cfa-45ed-b07e-78af2bbe3b35",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<html>\\(\\displaystyle \\left( D_{A}, D_{B} \\right) \\ {\\mapsto} \\ -\\frac{1}{200} \\, D_{A}^{2} - \\frac{1}{100} \\, {\\left(D_{A} - 10\\right)} D_{B} + \\frac{1}{200} \\, {\\left(D_{A} - 20\\right)} D_{B} + \\frac{1}{200} \\, D_{A} D_{B} + \\frac{1}{200} \\, D_{B}^{2} + \\frac{1}{10} \\, D_{A} - \\frac{1}{10} \\, D_{B} + \\frac{D_{B}}{2 \\, D_{A}} + \\frac{1}{2}\\)</html>"
],
"text/latex": [
"$\\displaystyle \\left( D_{A}, D_{B} \\right) \\ {\\mapsto} \\ -\\frac{1}{200} \\, D_{A}^{2} - \\frac{1}{100} \\, {\\left(D_{A} - 10\\right)} D_{B} + \\frac{1}{200} \\, {\\left(D_{A} - 20\\right)} D_{B} + \\frac{1}{200} \\, D_{A} D_{B} + \\frac{1}{200} \\, D_{B}^{2} + \\frac{1}{10} \\, D_{A} - \\frac{1}{10} \\, D_{B} + \\frac{D_{B}}{2 \\, D_{A}} + \\frac{1}{2}$"
],
"text/plain": [
"(D_A, D_B) |--> -1/200*D_A^2 - 1/100*(D_A - 10)*D_B + 1/200*(D_A - 20)*D_B + 1/200*D_A*D_B + 1/200*D_B^2 + 1/10*D_A - 1/10*D_B + 1/2*D_B/D_A + 1/2"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p(D_A,D_B) = integral(D_A*X_A/D_B*1/10, X_A, 0, D_B)/10 \\\n",
" + integral((10 - (D_A-D_B+X_A))/10, X_A, 0, D_B)/10 \\\n",
" + integral((D_A-D_B+X_A)/10, X_A, D_B, 10-(D_A-D_B))/10 \\\n",
" + integral(1, X_A, 10-(D_A-D_B), 10)/10 \\\n",
" + integral((10-D_A*X_A/D_B)/10, X_A, D_B, 10*D_B/D_A)/10\n",
"p"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "691d33bd-1e9c-433e-90d0-20e9dc75102a",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<html>\\(\\displaystyle \\left(\\frac{577}{600}, 0.961666666666667\\right)\\)</html>"
],
"text/latex": [
"$\\displaystyle \\left(\\frac{577}{600}, 0.961666666666667\\right)$"
],
"text/plain": [
"(577/600, 0.961666666666667)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(p(6,5), numerical_approx(p(6,5)))"
]
},
{
"cell_type": "markdown",
"id": "b567242b-a2bc-40d7-91a4-b392cea84230",
"metadata": {},
"source": [
"This is excellent agreement with the results of the simulations."
]
},
{
"cell_type": "markdown",
"id": "20bc4e52-ee4e-4cc3-9aba-4990c7500ce4",
"metadata": {},
"source": [
"Extra credit\n",
"------------\n",
"Integrating over $D_A$ and $D_B$, assuming $D_A > D_B$, then doubling to count the $D_A < D_B$ cases,"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "faef0c26-d79a-412d-b3d2-405b9a8c41aa",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<html>\\(\\displaystyle \\left(\\frac{11}{12}, 0.916666666666667\\right)\\)</html>"
],
"text/latex": [
"$\\displaystyle \\left(\\frac{11}{12}, 0.916666666666667\\right)$"
],
"text/plain": [
"(11/12, 0.916666666666667)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"extra_credit = 2*integral(integral(p, D_B, 0, D_A)/10, D_A, 0, 10)/10\n",
"(extra_credit, numerical_approx(extra_credit))"
]
},
{
"cell_type": "markdown",
"id": "84e16fcc-ee38-4b15-8565-5362b3fc203a",
"metadata": {},
"source": [
"This is also excellent agreement with the results of the simulations."
]
},
{
"cell_type": "markdown",
"id": "5b398e83-dda0-42bc-a4a7-fe4f7c94bc8e",
"metadata": {},
"source": [
"Making the rounds\n",
"-----------------\n",
"If an $n$-dimensional ball of radius $r$ has volume $A_nr^n$, then an $n+1$-dimensional ball will be made up\n",
"of $n$-dimensional balls of radius $\\rho_n$ that are centered $x_n$ from the origin along the axis in the\n",
"$n+1$th dimension, where $\\rho_n^2 + x_n^2 = 1$, so $\\rho_n = \\sqrt{1 - x_n^2}$. The volume is then\n",
"\n",
"$$\n",
" \\int_{-1}^1 dx\\, A_n\\rho^n = 2\\int_0^1 dx\\, A_n (1-x^2)^{n/2}\n",
"$$\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "fa5481f3-86e2-4670-9c64-2b976df56742",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<html>\\(\\displaystyle n \\ {\\mapsto}\\ \\operatorname{B}\\left(\\frac{1}{2}, \\frac{1}{2} \\, n + 1\\right)\\)</html>"
],
"text/latex": [
"$\\displaystyle n \\ {\\mapsto}\\ \\operatorname{B}\\left(\\frac{1}{2}, \\frac{1}{2} \\, n + 1\\right)$"
],
"text/plain": [
"n |--> beta(1/2, 1/2*n + 1)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = var('n')\n",
"assume(n > 0)\n",
"f(n) = 2*integral((1-x^2)^(n/2), x, 0, 1)\n",
"f"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "549707bf-69fd-4eea-aa5f-77a3bcb0a8da",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<html>\\(\\displaystyle n \\ {\\mapsto}\\ {\\prod_{j_{1}=1}^{n} \\operatorname{B}\\left(\\frac{1}{2}, \\frac{1}{2} \\, j_{1} + \\frac{1}{2}\\right)}\\)</html>"
],
"text/latex": [
"$\\displaystyle n \\ {\\mapsto}\\ {\\prod_{j_{1}=1}^{n} \\operatorname{B}\\left(\\frac{1}{2}, \\frac{1}{2} \\, j_{1} + \\frac{1}{2}\\right)}$"
],
"text/plain": [
"n |--> product(beta(1/2, 1/2*j1 + 1/2), j1, 1, n)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m = var('m')\n",
"A(n) = product(f(m), m, 0, n-1)\n",
"A"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "8a0604d5-4336-40c7-9882-1cfe59182d70",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<html>\\(\\displaystyle \\left[1, 2 \\, r, \\pi r^{2}, \\frac{4}{3} \\, \\pi r^{3}, \\frac{1}{2} \\, \\pi^{2} r^{4}, \\frac{8}{15} \\, \\pi^{2} r^{5}, \\frac{1}{6} \\, \\pi^{3} r^{6}, \\frac{16}{105} \\, \\pi^{3} r^{7}, \\frac{1}{24} \\, \\pi^{4} r^{8}, \\frac{32}{945} \\, \\pi^{4} r^{9}, \\frac{1}{120} \\, \\pi^{5} r^{10}, \\frac{64}{10395} \\, \\pi^{5} r^{11}, \\frac{1}{720} \\, \\pi^{6} r^{12}, \\frac{128}{135135} \\, \\pi^{6} r^{13}, \\frac{1}{5040} \\, \\pi^{7} r^{14}, \\frac{256}{2027025} \\, \\pi^{7} r^{15}\\right]\\)</html>"
],
"text/latex": [
"$\\displaystyle \\left[1, 2 \\, r, \\pi r^{2}, \\frac{4}{3} \\, \\pi r^{3}, \\frac{1}{2} \\, \\pi^{2} r^{4}, \\frac{8}{15} \\, \\pi^{2} r^{5}, \\frac{1}{6} \\, \\pi^{3} r^{6}, \\frac{16}{105} \\, \\pi^{3} r^{7}, \\frac{1}{24} \\, \\pi^{4} r^{8}, \\frac{32}{945} \\, \\pi^{4} r^{9}, \\frac{1}{120} \\, \\pi^{5} r^{10}, \\frac{64}{10395} \\, \\pi^{5} r^{11}, \\frac{1}{720} \\, \\pi^{6} r^{12}, \\frac{128}{135135} \\, \\pi^{6} r^{13}, \\frac{1}{5040} \\, \\pi^{7} r^{14}, \\frac{256}{2027025} \\, \\pi^{7} r^{15}\\right]$"
],
"text/plain": [
"[1,\n",
" 2*r,\n",
" pi*r^2,\n",
" 4/3*pi*r^3,\n",
" 1/2*pi^2*r^4,\n",
" 8/15*pi^2*r^5,\n",
" 1/6*pi^3*r^6,\n",
" 16/105*pi^3*r^7,\n",
" 1/24*pi^4*r^8,\n",
" 32/945*pi^4*r^9,\n",
" 1/120*pi^5*r^10,\n",
" 64/10395*pi^5*r^11,\n",
" 1/720*pi^6*r^12,\n",
" 128/135135*pi^6*r^13,\n",
" 1/5040*pi^7*r^14,\n",
" 256/2027025*pi^7*r^15]"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"r = var('r')\n",
"[A(i).simplify()*r^i for i in [0..15]]"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "436cfb86-753e-4f07-a405-dc67bc4c73c4",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<html>\\(\\displaystyle \\left[1.00, 2.00 \\, r, 3.14 \\, r^{2}, 4.19 \\, r^{3}, 4.94 \\, r^{4}, 5.26 \\, r^{5}, 5.17 \\, r^{6}, 4.72 \\, r^{7}, 4.06 \\, r^{8}, 3.30 \\, r^{9}, 2.55 \\, r^{10}, 1.88 \\, r^{11}, 1.34 \\, r^{12}, 0.911 \\, r^{13}, 0.599 \\, r^{14}, 0.381 \\, r^{15}\\right]\\)</html>"
],
"text/latex": [
"$\\displaystyle \\left[1.00, 2.00 \\, r, 3.14 \\, r^{2}, 4.19 \\, r^{3}, 4.94 \\, r^{4}, 5.26 \\, r^{5}, 5.17 \\, r^{6}, 4.72 \\, r^{7}, 4.06 \\, r^{8}, 3.30 \\, r^{9}, 2.55 \\, r^{10}, 1.88 \\, r^{11}, 1.34 \\, r^{12}, 0.911 \\, r^{13}, 0.599 \\, r^{14}, 0.381 \\, r^{15}\\right]$"
],
"text/plain": [
"[1.00,\n",
" 2.00*r,\n",
" 3.14*r^2,\n",
" 4.19*r^3,\n",
" 4.94*r^4,\n",
" 5.26*r^5,\n",
" 5.17*r^6,\n",
" 4.72*r^7,\n",
" 4.06*r^8,\n",
" 3.30*r^9,\n",
" 2.55*r^10,\n",
" 1.88*r^11,\n",
" 1.34*r^12,\n",
" 0.911*r^13,\n",
" 0.599*r^14,\n",
" 0.381*r^15]"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[numerical_approx(A(i).simplify(),digits=3)*r^i for i in [0..15]]"
]
},
{
"cell_type": "markdown",
"id": "4a01b7c0-07e0-4171-b3a4-f698ff320039",
"metadata": {},
"source": [
"Since $A_{n+1} = A_n \\text{B}\\left(\\frac{1}{2},\\frac{n+2}{2}\\right)$, and $\\text{B}(x,y) \\sim \\Gamma(x) y^{-x}$\n",
"for large $y$, then\n",
"\n",
"$$\n",
" A_{n+1} \\sim A_n \\Gamma\\left(\\frac12\\right) \\sqrt{\\frac{2}{n+2}} = A_n\\sqrt{\\frac{2\\pi}{n+2}}\n",
"$$\n",
"\n",
"so the volume of a unit ball goes to zero for large $n$, and the maximum is $8\\pi^2/15$ at $n = 5$, which\n",
"makes sense because $\\frac{2\\pi}{n+2} > 1$ for $n < 5$ and $\\frac{2\\pi}{n+2} < 1$ for $n \\ge 5$."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "a239a532-c122-43f9-8a13-41352d2243ec",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<html>\\(\\displaystyle \\left[\\left(1.88 \\, r^{11}, 1.85 \\, r^{11}\\right), \\left(1.34 \\, r^{12}, 1.31 \\, r^{12}\\right), \\left(0.911 \\, r^{13}, 0.895 \\, r^{13}\\right), \\left(0.599 \\, r^{14}, 0.589 \\, r^{14}\\right), \\left(0.381 \\, r^{15}, 0.376 \\, r^{15}\\right), \\left(0.235 \\, r^{16}, 0.232 \\, r^{16}\\right), \\left(0.141 \\, r^{17}, 0.139 \\, r^{17}\\right), \\left(0.0822 \\, r^{18}, 0.0811 \\, r^{18}\\right), \\left(0.0466 \\, r^{19}, 0.0460 \\, r^{19}\\right), \\left(0.0258 \\, r^{20}, 0.0255 \\, r^{20}\\right), \\left(0.0139 \\, r^{21}, 0.0138 \\, r^{21}\\right)\\right]\\)</html>"
],
"text/latex": [
"$\\displaystyle \\left[\\left(1.88 \\, r^{11}, 1.85 \\, r^{11}\\right), \\left(1.34 \\, r^{12}, 1.31 \\, r^{12}\\right), \\left(0.911 \\, r^{13}, 0.895 \\, r^{13}\\right), \\left(0.599 \\, r^{14}, 0.589 \\, r^{14}\\right), \\left(0.381 \\, r^{15}, 0.376 \\, r^{15}\\right), \\left(0.235 \\, r^{16}, 0.232 \\, r^{16}\\right), \\left(0.141 \\, r^{17}, 0.139 \\, r^{17}\\right), \\left(0.0822 \\, r^{18}, 0.0811 \\, r^{18}\\right), \\left(0.0466 \\, r^{19}, 0.0460 \\, r^{19}\\right), \\left(0.0258 \\, r^{20}, 0.0255 \\, r^{20}\\right), \\left(0.0139 \\, r^{21}, 0.0138 \\, r^{21}\\right)\\right]$"
],
"text/plain": [
"[(1.88*r^11, 1.85*r^11),\n",
" (1.34*r^12, 1.31*r^12),\n",
" (0.911*r^13, 0.895*r^13),\n",
" (0.599*r^14, 0.589*r^14),\n",
" (0.381*r^15, 0.376*r^15),\n",
" (0.235*r^16, 0.232*r^16),\n",
" (0.141*r^17, 0.139*r^17),\n",
" (0.0822*r^18, 0.0811*r^18),\n",
" (0.0466*r^19, 0.0460*r^19),\n",
" (0.0258*r^20, 0.0255*r^20),\n",
" (0.0139*r^21, 0.0138*r^21)]"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[(numerical_approx(A(i+1).simplify(),digits=3)*r^(i+1),numerical_approx(A(i).simplify()*sqrt(2*pi/(i+2)),digits=3)*r^(i+1)) for i in [10..20]]"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "SageMath 10.3",
"language": "sage",
"name": "sagemath-10.3"
},
"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.11.8"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment