{ "cells": [ { "metadata": {}, "cell_type": "markdown", "source": [ "# Hubble (1929): Expansion of the Universe\n", "\n", "Edwin Hubble is famous for his observation that the universe is expanding (or more technically that galaxies are moving away from Earth at speeds proportional to their distance). From this observation comes Hubble's Law and Hubble's Constant. In the example, we will repeat Hubble's analysis but using modern estimating equations\n", "\n", "The corresponding data for this analysis can be obtained from the R package `gamair`.\n", "\n", "## Setup" ], "id": "f32e93be52118619" }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-10T15:04:54.441851800Z", "start_time": "2026-04-10T15:04:54.336889500Z" } }, "cell_type": "code", "source": [ "import numpy as np\n", "import scipy as sp\n", "import pandas as pd\n", "import matplotlib\n", "import matplotlib.pyplot as plt\n", "\n", "import delicatessen as deli\n", "from delicatessen import MEstimator\n", "from delicatessen.estimating_equations import ee_regression, ee_robust_regression\n", "\n", "print(\"NumPy version: \", np.__version__)\n", "print(\"SciPy version: \", sp.__version__)\n", "print(\"Pandas version: \", pd.__version__)\n", "print(\"Matplotlib version: \", matplotlib.__version__)\n", "print(\"Delicatessen version:\", deli.__version__)" ], "id": "9028df9470c94e7a", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "NumPy version: 2.3.5\n", "SciPy version: 1.16.3\n", "Pandas version: 2.3.3\n", "Matplotlib version: 3.10.8\n", "Delicatessen version: 4.2\n" ] } ], "execution_count": 12 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-10T14:56:34.318425300Z", "start_time": "2026-04-10T14:56:34.293249600Z" } }, "cell_type": "code", "source": "d = pd.read_csv(\"data/hubble.csv\")", "id": "25b03e8a50abe06d", "outputs": [], "execution_count": 4 }, { "metadata": {}, "cell_type": "markdown", "source": [ "This data set is consists of observations of the distances (`x`) and velocities (`y`) for 24 galaxies containing Cepheid stars.\n", "\n", "The following code plots the observed data points" ], "id": "78d44365636033a0" }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-10T14:58:23.390767300Z", "start_time": "2026-04-10T14:58:23.161944200Z" } }, "cell_type": "code", "source": [ "plt.plot(d['x'], d['y'], 'o')\n", "plt.xlabel(\"Distance\")\n", "plt.ylabel(\"Velocity\")\n", "plt.tight_layout()" ], "id": "9a3a7623799b1175", "outputs": [ { "data": { "text/plain": [ "
" ], "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAHWCAYAAAARl3+JAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPk9JREFUeJzt3Qd8VFXe//FfQhLUmAAKJJE/IEVFIVQbKEUBy2MJ9roiuj6iPqsUhWQVQV0NuCtBQX0UEVFE15W6gjTlAUtEKUJAsNAXUmghESIt5//6HXdmM2mkTDJ37nzer9d9Ze69J5PbJvnmnHvODRMRIwAAAAhq4YHeAAAAAFQfoQ4AAMAFCHUAAAAuQKgDAABwAUIdAACACxDqAAAAXIBQBwAA4AKEOgAAABeICPQGBIszzjhD8vPzA70ZAADA5WJiYmTXrl2V/j5CXQUD3c6dO6tyXgAAACqtSZMmlQ52hLoK8NTQ6QGmtg4AANRkLZ1WJFUlbxDqKkEPMKEOAAA4ER0lAAAAXIBQBwAA4AKEOgAAABcg1AEAALgAoQ4AAMAFCHUAAAAuQKgDAABwAUIdAACACxDqAAAAXIBQBwAA4AI8JgwAEDLCwsOlZecOEtuooeTt3iObV60RU1gY6M0C/IJQBwAICYm9e0q/5MFSPz7Ouyw3K1tmjU6TjM+WBnTbAH+g+RUAEBKBrv/YVKnXuJHPcp3X5boeCHaEOgCA65tctYZOxNjXxdfp8qThg0qsA4INVzAAwNX0Hjptci0rtOnyBgnxthwQzAh1AABX004R/iwHOBWhDgDgatrL1Z/lAKci1AEAXE2HLdFermUNXaLL92dm2XJAMCPUAQBcTUObDlsiElYi2P0+Hyazx4xjvDoEPUIdAMD1dBy6KUNS5EDObp/ludk5djnj1MEtTKCm7t27mzlz5pidO3calZSU5LO+LI8//ri3zJYtW0qsHz58uM/7JCYmmmXLlpmCggKzfft288QTT1RqO2NiYuz76tdAHi8mjgHXANcA10D1roGw8HDT6vxOptPVfe1XneeY8rkSBx2D6mSOgD5RIjo6WtasWSNvv/22zJw5s8T6+Ph4n/mrr75aJk2aJNOnT/dZPmLECJk4caJ3Pj8/3/s6JiZGFi5cKIsXL5aBAwdKYmKi/Xm5ubk+3wMAcD9tbt20YnWgNwOoEQENdfPnz7dTWbKzs33mk5KSZMmSJbJlyxaf5Rriipf1uOuuuyQqKkruu+8+OXr0qPzwww/SsWNHGTJkCKEOAAC4RtDcU9e4cWO55pprbE1dccnJybJnzx5ZtWqVPP7441KnTh3vuq5du8qyZctsoPNYsGCBtGnTRurXr19r2w8AAODamrrK6N+/v62RmzFjhs/yV155xYa5ffv2Sbdu3SQ1NVUSEhJk6NCh3ibc4jV7nlo9XafNsMVpzV7dunV9mnABAACcLGhCnTafvv/++3L48GGf5Wlp2k39dxkZGXLkyBF54403JCUlxb6uCv3eUaNGVXubAQAAaktQNL9eeumltrn0rbfeOmHZ5cuXS2RkpJx55pl2PisrS+Li4nzKeOZ1XWm0ti82NtY7NWnSxC/7AQAAENKh7v7775cVK1bI2rVrT1hWO0EcP35ccnJy7Hx6err06NFDIiL+UynZt29f2bhxY6lNr0pr+LSpt+gEAADgdAEbiyU6Otp06NDBTmrQoEH2ddOmTX3Ga/n111/Ngw8+WOL7L774YvPYY4+Z9u3bmxYtWpg777zTZGdnm3feecdbJjY21mRmZpopU6aY8847z9x66632/R544IFaGTOGiWPANcA1wDXANcA1wDUgtZM5Anege/bsWergwpMnT/aW0fB18OBBG86Kf3+nTp1Menq62b9/vzl06JBZv369SU5ONlFRUWUOPrxjxw4zbNiwSm0noY4PI7+QuQa4BrgGuAa4BqQWjkF1MkeYJ9mhbNr7NS8vz95fR1MsAABwYuYIinvqAAAAUD5CHQAAgAsQ6gAAAFyAUAcAAOAChDoAAAAXINQBAAC4AKEOAADABQh1AAAALkCoAwAAcAFCHQAAgAsQ6gAAAFyAUAcAAOAChDoAAAAXINQBAAC4AKEOAADABQh1AAAALkCoAwAAcAFCHQAAgAsQ6gAAAFyAUAcAAOAChDoAAAAXINQBAAC4AKEOAADABSICvQEAAKBsYeHh0rJzB4lt1FDydu+RzavWiCks5JChBEIdAAAOldi7p/RLHiz14+O8y3KzsmXW6DTJ+GxpQLcNzkPzKwAADg10/cemSr3GjXyW67wu1/VAUYQ6AAAc2OSqNXQixr4uvk6XJw0fVGIdQhtXAwAADqP30GmTa1mhTZc3SIi35QAPQh0AAA6jnSL8WQ6hgVAHAIDDaC9Xf5ZDaCDUAQDgMDpsifZyLWvoEl2+PzPLlgM8CHUAADiMhjYdtkQkrESw+30+TGaPGcd4dfBBqAMAwIF0HLopQ1LkQM5un+W52Tl2OePUobgw2y8a5YqJiZG8vDyJjY2V/Px8jhYAoNbwRInQElONzMETJQAAcDBtbt20YnWgNwNBgOZXAAAAFyDUAQAAuAChDgAAwAUIdQAAAC5AqAMAAHABQh0AAIALEOoAAABcgFAHAADgAgENdd27d5c5c+bIzp07xRgjSUlJPusnT55slxedPv30U58yDRo0kKlTp8qBAwdk//798tZbb0l0dLRPmcTERFm2bJkUFBTI9u3b5YknnqiV/QMAAAiJUKfha82aNfLII4+UWUZDXHx8vHe64447fNa///770rZtW+nbt69ce+210qNHD3nzzTd9HrexcOFC2bZtm3Tp0sUGulGjRskDDzxQo/sGAABQ24wTJpWUlOSzbPLkyWbmzJllfk+bNm3s93Xp0sW77MorrzTHjx83CQkJdn7gwIFm7969JjIy0lsmNTXVbNiwocLbFhMTY3+Ofg30cWLiGHANcA1wDXANcA249xqIqUbmcPw9db169ZLs7GzZuHGjvPbaa3Laaad513Xt2tU2ua5cudK7bPHixVJYWCgXXXSRt4w2vR49etRbZsGCBdKmTRupX79+Le8NAABAzYiooff1i/nz58uMGTNky5Yt0qpVK3nhhRdsc6wGNQ1u2hybk5Pj8z3Hjx+Xffv22XVKv+r3F6Uh0bMuNze3xM+NioqSunXr+jThAgAAOJmjQ93f//537+t169bJ2rVrZfPmzbb27vPPP6+xn5uSkmLvuwMAAAgWjm9+LUpr3Hbv3i2tW7e281lZWdK4cWOfMnXq1LFNtLrOUyYuLs6njGfeU6a41NRUiY2N9U5NmjSpoT0CAABOFhYeLq3O7ySdru5rv+q8Uzm6pq44DVenn366ZGZm2vn09HQ7pEnnzp1l1apVdtnll18u4eHhsnz5cm+Z559/XiIiIuTYsWN2mfaU1Xv0Smt6VUeOHLETAAAIXYm9e0q/5MFSP/4/lUO5Wdkya3SaZHy2VJwm4EOadOjQwU6qRYsW9nXTpk3tuhdffNF2eGjevLkNa7Nnz5ZffvnFdnRQGsz0HruJEyfKBRdcIN26dZMJEybIhx9+6A1+06ZNswFt0qRJct5558mtt94qjz32mIwdOzaQuw4AABwe6PqPTZV6jRv5LNd5Xa7rnShg3XZ79uxpSqNDmZx00klm/vz5Jjs72xw+fNhs2bLFvPHGG6Zx48Y+79GgQQPz/vvvm7y8PJObm2smTZpkoqOjfcokJiaaZcuWmYKCArNjxw4zbNiwWutezMQx4BrgGuAa4BrgGgiuayAsPNyMWDTL/G3NV+aljPQSky5/auFMW87fP7s6mSPs3y9QDu39mpeXZ++vy8/P51gBAOBirc7vJA9Pfu2E5V4b8LBsWrHaMZnDuXf7AQAABEBso4Z+LVdbgqqjBAAAZdFeiS07d7B/aPN275HNq9aIKSzkgKHS9PrxZ7naQqgDAAS9YOulCGfbvGqNvX60U0RpQ5joPwu52Tm2nJPQ/AoACOpxwYK1lyKcyxQW2n8IRMJK1Pb+Ph8ms8eMc1xNMB0lKoCOEgDgzBo3DXxPLZhxwhqV56+6yXF/gBGc1+P+zCwb6GqqBrg6mYNQV8MHGABQsXvfPDVuOihD0YDmqRmZMiSlxB/SQPZSRGgIq+V7NauTObinDgAQ0JqPX/ftk+nPvyRJTzxaItApndc/oknDB8m6JV/4/EEN1l6KCB6msDBo/iHgnjoAQK0o6963U087Te75219s0Cvr/jld3iAh3taYuKGXIlATCHUAgBqnoUxr6EqriauM4jVunl6KZTWH6XK9B8ppvRSBmkCoAwDUOK1hK7cmLkxv8a58jVuw9lIEagKhDgBQ46p7T1t5NW7aeUI7URzI2e2zXHu9lta5Aqit4XZqGx0lAAA1rjL3tGmAK633a3k1bhrctBMFT5RAKA9wzZAmFcCQJgBQzT824eEyask/baeIE9HesEXL1fS4YMCJVGW4napiSBMAgKPpH7/pz/1V7hn7Qpn30BljbM3HC/91i7TomMgzXOH4Tj5h5Qy3Ewg0vwIAasXaxf8nSyZPlcsG3F1qoBMjtkau8NixGh8XrLYHlEXwd/IpS9HhdgI9nh2hDgBQa+amvSbb122Qm596Qk49rYF3udbQ1VYTazDcGwXniA2iAa4JdQCAWpWxaIms+2xpQGrKit4bVZQOiKzL6S2L4oJpgGtCHQAgJB69FEz3RsE5Nv97gGsN/qUNYaLXig6f44QBrp01wAoAAIEaALmMR5EhtJkgGuCaUAcACAnBdG8UnCUjSAa4pvkVABASguneKDhPRhAMcE2oAwCEhGC6NwrOZAJwL2hl0PwKAAgJwXRvFFAVhDoAQMgIlnujgKrg2a8VwLNfAcBdeKIEnIpnvwIA4KJ7o4CqoPkVAADABQh1AAAALkCoAwAAcAFCHQAAgAsQ6gAAAFyAUAcAAOAChDoAAAAXINQBAAC4AKEOAADABQh1AAAALkCoAwAAcAFCHQAAgAsQ6gAAAFyAUAcAAOAChDoAAAAXINQBAAC4AKEOAADABQIa6rp37y5z5syRnTt3ijFGkpKSvOsiIiJk9OjRsnbtWvn1119tmSlTpkhCQoLPe2zZssV+b9Fp+PDhPmUSExNl2bJlUlBQINu3b5cnnnii1vYRAADA9aEuOjpa1qxZI4888kiJdaeccop07txZnnvuOfv1xhtvlHPOOceGwOJGjBgh8fHx3mn8+PHedTExMbJw4ULZtm2bdOnSxQa6UaNGyQMPPFDj+wcAAFCbjBMmlZSUVG6Z888/35Zr2rSpd9mWLVvMY489Vub3DBw40Ozdu9dERkZ6l6WmppoNGzZUeNtiYmLsz9WvgT5OTBwDrgGuAa4BrgGuAfdeAzHVyBxBdU9dvXr1pLCwUHJzc32WJycny549e2TVqlXy+OOPS506dbzrunbtaptejx496l22YMECadOmjdSvX79Wtx8AAKCmREiQqFu3rowZM0Y++OADyc/P9y5/5ZVXbJjbt2+fdOvWTVJTU+19d0OHDrXrtTlW77srKjs727uueEBUUVFR9ucVbcIFAABwsqAIddpp4qOPPpKwsDB56KGHfNalpaV5X2dkZMiRI0fkjTfekJSUFPu6KvR79b47AACAYBEeLIGuefPm0rdvX59autIsX75cIiMj5cwzz7TzWVlZEhcX51PGM6/rSqO1fbGxsd6pSZMmftsfAACAkAt1nkB31llnSZ8+fWwT64l07NhRjh8/Ljk5OXY+PT1devToYd/LQ8Phxo0bS216VVrDp+Gx6AQAAOBkEYEe0qR169be+RYtWkiHDh1seMvMzJSPP/7YDmdy7bXX2s4Pnho2Xa8dHy6++GK56KKLZMmSJTZ4aacIbY6dOnWqN7BNmzZNRo4cKZMmTbL35LVr104ee+wxGTx4cMD2GwAAoCYErNtuz549TWkmT55smjdvbsqi36ff36lTJ5Oenm72799vDh06ZNavX2+Sk5NNVFSUz89JTEw0y5YtMwUFBWbHjh1m2LBhtda9mIljwDXANcA1wDXANcA1ILWQOcL+/QLl0N6veXl59v46mmIBAIATM4ej76kDAABAxRDqAAAAXIBQBwAA4AKEOgAAABcg1AEAALgAoQ4AAMAFCHUAAAAuQKgDAABwAUIdAACACxDqAAAAXIBQBwAA4AKEOgAAABcg1AEAALgAoQ4AAMAFCHUAAAAuQKgDAABwAUIdAACACxDqAAAAXIBQBwAA4AKEOgAAABcg1AEAALhARKA3AABCVVh4uLTs3EFiGzWUvN17ZPOqNWIKCwO9WQCCFKEOAAIgsXdP6Zc8WOrHx3mX5WZly6zRaZLx2VLOCYBKo/kVAAIQ6PqPTZV6jRv5LNd5Xa7rAaCyCHUAUMtNrlpDJ2Ls6+LrdHnS8EEl1gHAifBbAwBqkd5Dp02uZYU2Xd4gId6WA4DKINQBQC3SThH+LAcAHoQ6AKhF2svVn+UAwINQBwC1SIct0V6uZQ1dosv3Z2bZcgBQGYQ6AKhFGtp02BKRsBLB7vf5MJk9Zhzj1QGoNEIdANQyHYduypAUOZCz22d5bnaOXc44dQCqIsz2n0e5YmJiJC8vT2JjYyU/P5+jBcAveKIEAH9mDp4oAQABos2tm1as5vgD8AuaXwEAAFyAUAcAAOAChDoAAAAXINQBAACEYqjbsmWLjBgxQpo2bVozWwQAAICaD3Xjxo2TG2+8UTZv3iwLFy6U2267TaKioir/kwEAABC4UPfyyy9Lp06d5MILL5QNGzbI+PHjJTMz037V5QAAAAgMU50pIiLCPProo6agoMAcO3bMrF692gwYMKBa7+m0KSYmxij9GuhtYeIYcA1wDXANcA1wDbj3GoipRuao8uDDERERcsMNN8iAAQOkb9++8s0338ikSZPk//2//ycvvPCC9OnTR+666y7/xk8AAACUqtKhTptYNcjdcccdUlhYKO+++64MHjxYfvzxR2+ZmTNnynfffVfZtwYAAEBt3VOnYe2ss86Shx56SJo0aSJPPPGET6Dz9JD98MMPT/he3bt3lzlz5sjOnTvFGCNJSUklyjzzzDOya9cuOXTokCxatEhat27ts75BgwYydepUOXDggOzfv1/eeustiY6O9imTmJgoy5Ytk4KCAtm+fbvdZgBA7T3jttX5naTT1X3tV50HUDMq1V7brFkzv7UbX3XVVea5554z/fr1s+3HSUlJPuuHDRtm9u/fb66//nqTmJhoZs2aZTZt2mTq1q3rLTNv3jx7H9+FF15oLrnkEvPTTz+Z999/36dtOjMz07z33nvmvPPOM7fddps5ePCgeeCBB2qlfZuJY8A1wDUQytdAYu+eZsSiWealjHTvpPO6PNDbxsQxEAceg2pmjsp9g4aq0047rcTyevXq2XVV3YnSQt2uXbvM0KFDvfOxsbG2Q4YGM51v06aN/b4uXbp4y1x55ZXm+PHjJiEhwc4PHDjQ7N2710RGRnrLpKammg0bNtTWAWbiGHANcA2E5DWgwe1va742f1vzlU+o03ldTrAL/DliEscdg+pkjkrXgZ955plSp06dEsvr1q1rm2P9pUWLFpKQkCCLFy/2LsvLy5Ply5dL165d7bx+1SbXlStXestoeb3X76KLLvKW0abXo0ePesssWLBA2rRpI/Xr1/fb9gIA/kObWPslD7Z/a4o3t/4+byRp+CCaYoFAdJS47rrrvK+vvPJKew+bh4a83r17y9atW/22YfHx8fZrdna2z3Kd96zTrzk5OT7rjx8/Lvv27fMpo/f4FX8Pz7rc3NwSP1sHU9aQ6hETE+O3/QKAUNCycwepHx9X5noNdg0S4m25TStW1+q2ARLqoW7WrFn2q3ZomDJlis86rQXTQDd06FBxg5SUFBk1alSgNwMAglZso4Z+LQfgxCrc/Kq1cTpp79HGjRt753U66aSTbHPm3LlzxV+ysrLs17g43//0dN6zTr/qthTfztNOO82nTGnvUfRnFJeamiqxsbHeyZ/NygAQCvJ27/FrOQAnVul76lq2bCl79+6VmqZNpvr4MW3WLdoMqvfKpaen23n9qkOadO7c2Vvm8ssvl/DwcHvvnadMjx497GDJHjpY8saNG0ttelVHjhyR/Px8nwkAUHGbV62R3KxsMYWFpa7X5fszs2w5ALXY/PqnP/1J3nzzTTl8+LB9XR59BmxF6XhyRced084RHTp0sPfE7dixQ8aNGydPPfWU/PzzzzbkPffcc3bMOk9TsAazTz/9VCZOnCgDBw6UyMhImTBhgh0jTwOhmjZtmowcOdI+7WLMmDHSrl07eeyxx+yAyQCAmqGhbdboNOk/NtW+LtpZ4vegFyazx4wrM/QBqJoTdpHdvHmzdxgTfV3WVNkhTXr27GlKM3nyZG+ZZ555xo4zp0OZLFq0yJx11lk+79GgQQM7Ll1eXp7Jzc01kyZNMtHR0b7d6hMTzbJly+x77Nixw45/V1vdi5k4BlwDXAOhfA2UNk7dUwtnMpyJA84NkzjyGFQnc4T9+wXKoc2+OpyK3l9HUywAVI7W0mkvV+0UoffQaZMrNXSA/zNHpZ/9CgBAZWiAY9gSwIEdJT7++GMZNmxYieX6PNWPPvrIX9sFAACAmgx12pN03rx5JZZrhwVdBwAAgCAIdaeeeqod8qM4HYBY238BAAAQBKEuIyNDbrvtthLLb7/9dvnhhx/8tV0AAACoyY4SOlbcjBkzpFWrVvL555/bZTpA8B133CG33HJLZd8OAAAAgQh1n3zyifTr10/+/Oc/y8033ywFBQWydu1a6dOnjyxbtswf2wQAAIBKYpy6CmCcOgAA4Npx6vR5q+eee659vX79evn++++r+lYAAACopkqHukaNGtlnq/bq1Utyc3Ptsvr168uSJUtsZ4k9e/ZUd5sAAABQ071fx48fb6sG27ZtK6effrqd2rVrZ6sJX3nllcq+HQAAAAJxT53WzmmniBUrVvgsv+CCC2ThwoXSoEEDcRvuqQMAAE7PHJWuqQsPD7cDDReny3QdAAAAal+lU5iOTffyyy9LQkKCd9kZZ5whaWlp8tlnn/l7+wAAAFATHSX+53/+R+bMmSNbt26VHTt22GVNmzaVdevWyd13313ZtwOAoBYWHi4tO3eQ2EYNJW/3Htm8ao2YwsJAbxaAEFTpUPevf/3LDmei99W1adPGLtuwYQO1dABCTmLvntIvebDUj4/zLsvNypZZo9Mk47OlAd02AKGHwYcrgI4SAEoLdP3Hptq+Zlpb5/F7LV2YTBmSQrAD4LzBh//0pz9VasgTAHAzDXFaQ1c80HnWabBLGj5I1i35gqZYALWmQqFu8GD95XVixhhCHQDX03voija5FqfBrkFCvC23acXqWt02AKGrQqGuZcuWNb8lABAktFOEP8sBgD9UeWC5yMhIOfvss6VOnTp+2RAACBbay9Wf5QAgIKHu5JNPlrfeeksOHTok69evl2bNmtnl+oiw4cOH+2WjAMDJdNgS7eVa1tAlunx/ZpYtBwCODXWpqanSoUMH6dWrl/z222/e5YsXL5bbbrvN39sHAI6joU2HLdFersWDnaf36+wx4+gkAcDZoa5fv352AOKvvvrKdozw0Fq7Vq1a+Xv7AMCRdBw6HbbkQM5un+W52TkMZwIgOAYfbtSokeTk5JRYHh0d7RPyACAUgp0OW8ITJQAEZahbsWKFXHPNNTJhwgQ77wlyf/zjHyU9Pd3/WwgADqbNrQxbAiCoQl3btm1tE2tKSorMnz9fzjvvPNsD9rHHHrOvu3XrJj179qzZrQUAAED17qlbu3atfPPNNzbAXXLJJRIREWGXXXHFFbY5tmvXrrJq1aqKvh0AAAACUVOntXADBgyQl156ScLDw2X69Ony+OOPyxdffOHP7QEAAEBN1tR9+eWXcv/990tCQoJ9FuyZZ54p//d//yc//vijDBs2TOLiyn5kDgAAAGpWmH0idRXpECZae/eHP/xB4uPj7b12SUlJ4jYxMTGSl5cnsbGxkp+fH+jNAQAALhVTjcxRrVCnTjnlFLnrrrvsoMT169e399q5DaEOAAA4PXNUOYF1795d7rvvPrnpppuksLBQPvroI5k0aVJV3w4AAADVUKlQp/fT3XvvvXZq3bq1fP311/Loo4/aQKfPggUAAIDDQ928efOkT58+smfPHnn33Xfl7bfflp9++qlmtw4AAAD+DXVHjx6Vm2++WT755BPb3AoAAIAgDHVu7NUKAAAQcuPUAQAAwLncN/4IAKBKwsLDpWXnDhLbqKHk7d4jm1etEcPtNkDQINQBcAUCSfUk9u4p/ZIHS/34/zwdKDcrW2aNTpOMz5ZW+/wAqHnVHnw4FDD4MOBsBJLqH7/+Y1PtnwMNxx6/19KFyZQhKQQ7IAgyB/fUAXBFIKnXuJHPcp3X5boeZdMQpzV0xQOdZ50uTxo+qMQ6AM7DpxRA0CKQVJ/eQ6dNrmWFNl3eICHelgPgbI4PdVu2bBFjTIlpwoQJdv2SJUtKrHv99dd93qNp06Z2fL2DBw9Kdna2vPjii1KnTp0A7REAfyGQVJ92ivBnOQCB4/iOEhdccIFPAGvXrp0sXrxY/vGPf3iXvfnmm/L0009754s+siw8PFzmzp0rWVlZ0q1bN/uoM30ihg6m/OSTT9bingDwNwJJ9WkvV3+WAxA4jg91+liyopKTk+WXX36RpUuX+oQ4rYErzRVXXCHnnXeefcRZTk6OrFmzRkaMGCFjxoyRUaNG2XAHIDgRSKpPhy3RXq56D2JpTbDaWSI3O8eWA+Bsjm9+LSoyMlLuvvtu+9zZou666y7ZvXu3ZGRkyAsvvCAnn3yyd13Xrl3tcg10HgsWLJB69epJ27Zta3X7AdRMIClrLDVdvj8zi0BSDj1GOmyJ9nItfhw9vV9njxnHeHVAEAiqUNevXz+pX7++vPPOO95l06ZNs0Hvsssuk9TUVPnDH/4gU6dO9a6Pj48vUYvnmdd1pYmKirJdiotOAJyHQOIfOg6dDltyIGe3z3KtoWM4EyB4OL75taj7779fPv30U8nMzPQumzhxovf1unXr7LrPP/9cWrZsKZs3b67Sz0lJSbFNswCCJ5CUGDg3O8fWMDFwbsWP47olX/BECSDImWCYmjVrZo4dO2auv/76csudcsopRl1xxRV2/plnnjGrV6/2KXPmmWfaMh07diz1PaKiokxMTIx3OuOMM2x5fR3o48DEMeAaKP0aCAsPN63O72Q6Xd3XftV5jhWfF64BrgEJsmOgWaOqmSNoauoGDBhg74vTnqzl6dixo/3qqc1LT0+3vVwbNWpk77tTffv2lQMHDsgPP/xQ6nscOXLETgCCqyl204rV1XoPHjUGIJgFRagLCwuzoW7KlCly/Phx73JtYr3zzjtl3rx5snfvXmnfvr2kpaXZnrHaOUItXLjQhrf33ntPhg0bZu+j+8tf/iKvvvoqwQ2AF48aAxDsgqKjhA5H0rx58xK9XrU2TddpcNu4caO89NJLMn36dLnuuuu8ZQoLC+Xaa6+1YVBr7bQThY5TV3RcOwChjUeNAXCDsH+3w6KGHq4LwNm0yfWpBTNOOE7b81fdxLAeABydOYKipg4AagqPGgPgFoQ6ACGNR40BcAtCHYCQxqPGALgFoQ5ASONRYwDcglAHIKTxqDEAbkGoAxDyePYpADdgSJMKYEgTIDTwRAkAwZw5guKJEgCcy01ByB+PGgOAQCHUAagyHq0FAM7BPXUAqoRHawGAsxDqAFSpybVf8mD7lMHij9b6fd5I0vBBpT52CwBQM/iNC6DSeLQWADgP99QBcNyjtdzU+QIAaguhDoCjHq1F5wsAqBqaXwE45tFadL4AgKoj1AFwxKO16HwBANVDqAPgiEdr0fkCAKqHe+oAVJkGt3VLvvBLp4aa7nwBOB0dhFBdhDoAjni0Vk12vgCcjg5C8AeaXwE//pfd6vxO0unqvvYrA+86o/MF4HR0EIK/UFMH+AH/Zfuv80X/san2ddFQXNXOF4DTnaiDkF7v+nQWvc2Bax8nQk0dUE38l+3czheA09FBCP5ETR1QDfyX7ezOF4DT0UEI/kSoA/zwX3Z5oa9BQrwt54/OBKHCX50vAKejgxD8ieZXoBr4LxtAddBBCP5EqAOqgf+yATjt6SwIXYQ6oBr4LxtAddFBCP7CPXVANTAMBwB/oIMQ/CHMDo6DcsXExEheXp7ExsZKfn4+RwsVGqdOB8rVZhOG4QAA1EbmINTV8AFG6OC5jQCAQGYOml8BP2EYDgBAINFRAgAAwAUIdQAAAC5AqAMAAHABQh0AAIALEOoAAABcgFAHAADgAoQ6AAAAFyDUAQAAuAChDgAAwAUIdQAAAC5AqAMAAHABQh0AAIALODrUjRw5UowxPtOGDRu86+vWrSsTJkyQPXv2SH5+vnz88cfSuHFjn/do2rSpfPLJJ3Lw4EHJzs6WF198UerUqROAvQEAAKg5EeJw69atkz59+njnjx075n2dlpYm11xzjdxyyy1y4MABG/BmzJghl156qV0fHh4uc+fOlaysLOnWrZskJCTIu+++K0ePHpUnn3wyIPsDAABQU4xTp5EjR5rVq1eXui42NtYcPnzY3HTTTd5l55xzjlEXXXSRnb/qqqvMsWPHTOPGjb1lHnzwQZObm2siIyMrvB0xMTH2ffVroI8JE8eAa4BrgGuAa4BrwL3XQEw1Moejm1/VWWedJTt37pRNmzbJ1KlTbXOq6tKli0RFRcnixYu9ZX/88UfZtm2bdO3a1c7r14yMDMnJyfGWWbBggdSrV0/atm0bgL0BAAAIwebX5cuXy7333mvDmjad6j12X3zxhbRr107i4+Pl8OHDttm1KL1vTtcp/arzxdd71pVFw6Ler+cRExPj5z0DAAAIoVA3f/5872utcdOQpzVxt956qxQUFNTYz01JSZFRo0bV2PsDAAD4m+ObX4vSWrmffvpJWrdubTs/aG2aNqUWFRcXZ9cp/arzxdd71pUlNTVVYmNjvVOTJk1qZH8AAABCMtRFR0dLq1atJDMzU1auXClHjhyR3r17e9efffbZ0rx5c0lPT7fz+jUxMVEaNWrkLdO3b18bDn/44Ycyf46+rw6RUnQCAABwOuPU6a9//avp0aOHad68uenatatZuHChycnJMQ0bNrTrX3vtNbN161bTq1cv07lzZ/PVV1/ZydsLJDzcrF271syfP9+0b9/eXHHFFSY7O9s8//zztdYThYljwDXANcA1wDXANcA1ILWTOZx7oD/44AOzc+dO89tvv5kdO3bY+ZYtW3rX161b10yYMMHs3bvX/Prrr2b69OkmLi7O5z2aNWtm5s6daw4ePGgDoQbFOnXqEOoccH6ZOAZcA1wDXANcA1wD4rdQF+ZJdiib9n7Ny8uz99fRFAsAAJyYOYLqnjoAAACUjlAHAADgAoQ6AAAAFyDUAQAAuAChDgAAwAUIdQAAAC5AqAMAAHABQh0AAIALEOoAAABcgFAHAADgAoQ6AAAAFyDUAQAAuAChDgAAwAUIdQAAAC5AqAMAAHABQh0AAIALEOoAAABcgFAHAADgAoQ6AAAAFyDUAQAAuAChDgAAwAUIdQAAAC5AqAMAAHABQh0AAIALEOoAAABcgFAHAADgAoQ6AAAAFyDUAQAAuEBEoDcAoSksPFxadu4gsY0aSt7uPbJ51RoxhYWB3iwAAIIWoQ61LrF3T+mXPFjqx8d5l+VmZcus0WmS8dlSzggAAFVA8ytqPdD1H5sq9Ro38lmu87pc1wMAgMoj1KFWm1y1hk7E2NfF1+nypOGDSqwDAAAnxl9P1Bq9h06bXMsKbbq8QUK8LQcAACqHUIdao50i/FkOAAD8B6EOtUZ7ufqzHAAA+A9CHWqNDluivVzLGrpEl+/PzLLlAABA5RDqUGs0tOmwJSJhJYLd7/NhMnvMOMarAwCgCgh1IUw7JrQ6v5N0urqv/VobvU51HLopQ1LkQM5un+W52Tl2OePUAQBQNWF2HAmUKyYmRvLy8iQ2Nlby8/NdcbQCPQAwT5QAAMC/mYNQF4KhzjMAcPHx4jxNoNSYAQAQfJmD5tcQwwDAAAC4E6EuxDAAMAAA7kSoCzEMAAwAgDs5OtQlJyfLt99+a9uWs7OzZebMmXL22Wf7lFmyZIkYY3ym119/3adM06ZN5ZNPPpGDBw/a93nxxRelTp06EooYABgAAHdydKjr2bOnvPrqq3LxxRdL3759JTIyUhYuXCinnHKKT7k333xT4uPjvdOwYcO868LDw2Xu3LkSFRUl3bp1k/79+8u9994rzz77rIQiBgAGAMCdgqr3a8OGDWX37t3So0cP+eKLL7w1dd9//70MHjy41O+56qqrbC3dGWecITk5OXbZgw8+KGPGjJFGjRrJ0aNHT/hz6f0KAABqQ8j0fq1Xr579um/fPp/ld911lw17GRkZ8sILL8jJJ5/sXde1a1e73BPo1IIFC+x7tW3bVkLRuiVfyILXJsqhPN+LhQGAAQAIXhESJMLCwmTcuHHy5Zdfyvr1673Lp02bJtu2bZNdu3ZJ+/btbQ3cOeecIzfddJNdr82xeh9dUZ55XVcabaqtW7euT2p2i9IGHT6Ye0C+mPp3WTxxCo/oAgAgSAVNqNN769q1ayeXXnqpz/KJEyd6X69bt04yMzPl888/l5YtW8rmzZur9LNSUlJk1KhR4jZFBx0u6pTYGLny4Qck65fNPKYLAIAgFRTNr+PHj5drr71WLrvsMtm5c2e5ZZcvX26/tm7d2n7NysqSuLj/1Eopz7yuK01qaqpty/ZMTZo0kWDHoMMAALhbeDAEuhtuuEEuv/xy2bp16wnLd+zY0X7VGjuVnp4uiYmJtlOEh/akPXDggPzwww+lvseRI0fszYlFp2DHoMMAALhbhNObXO+8805JSkqywcpTw6aB7LfffrNNrLp+3rx5snfvXntPXVpamixdutR2jlA6BIqGt/fee88OdaL30f3lL3+x763hLVQw6DAAAO7m6Jq6hx9+WOrXr29DmjaVeqbbbrvNrtdQ1qdPHxvcNm7cKC+99JJMnz5drrvuOu97FBYW2qbb48eP21q7qVOnyrvvvitPP/20hBIGHQYAwN2Capy6QHHDOHV6T91TC2ZIvcaN7OviTGGhHdLk+atuogcsAAABEjLj1KHqNLTNGp1mc7y+Lr5Ol88eM45ABwBAkCLUhZCMz5bKlCEpciBnt89yBh0GACD40fwaIs2vRWnzq/aG1c4Teq+dPg+2eO0dAAAIrszh6N6vqBka4DatWM3hBQDARWh+BQAAcAFCHQAAgAsQ6gAAAFyAUAcAAOAChDoAAAAXINQBAAC4AKEOAADABQh1AAAALkCoAwAAcAFCHQAAgAvwmDAH4FmsAACgugh1AZbYu6f0Sx4s9ePjvMtys7Jl1ug0yfhsaUC3DQAABA+aXwMc6PqPTZV6jRv5LNd5Xa7rAQAAKoJQF8AmV62hEzH2dfF1ujxp+KAS6wAAAEpDYgiQlp072CbXskKbLm+QEG/LAQAAnAihLkBiGzX0azkAABDaCHUBkrd7j1/LAQCA0EaoC5DNq9bYXq6msLDU9bp8f2aWLQcAAHAihLoA0dCmw5aIhJUIdr/Ph8nsMePKDH0AAABFEeoCSMehmzIkRQ7k7PZZnpudY5czTh0AAKioMDt2BsoVExMjeXl5EhsbK/n5+X4/WjxRAgAAVDdz8EQJB9Am1k0rVotTEToBAHA+Qh3KxWPMAAAIDtxThzLxGDMAAIIHoQ6l4jFmAAAEF0IdSsVjzAAACC6EOpSKx5gBABBcCHUoFY8xAwAguBDqUCoeYwYAQHAh1KFUPMYMAIDgQqhDmXiMGQAAwYPHhDngMWFOxxMlAACoHTwmDCH9GDMAAEDzKwAAgCtwTx0AAIALEOoAAABcgFAHAADgAoQ6AAAAFyDUAQAAuEBIhbqHH35YtmzZIgUFBfLNN9/IBRdcEOhNAgAA8IuQCXW33nqrjB07Vp555hnp3LmzrFmzRhYsWCCNGjUK9KYBAABUW8iEuiFDhsjEiRPlnXfekQ0bNsjAgQPl0KFDct999wV60wAAAKotJEJdZGSkdOnSRRYvXuxdZoyx8127di1RPioqyj6mo+gEAADgZCER6ho2bCgRERGSnZ3ts1zn4+PjS5RPSUmxz3r1TDt37qzFrQUAAKi8iCp8j+ulpqba++88tKZOgx01dgAAoCZVJ2uERKjbs2ePHDt2TOLi4nyW63xWVlaJ8keOHLFT8QNMjR0AAKgNmj3y8/Mr9T0hEeqOHj0qK1eulN69e8vs2bPtsrCwMDs/YcKEE37/rl27pEmTJpU+uDXJU3votO2qaaG436G4z4r9Dp3zzbnmXLtdTCV/j2t5zR6VFRKhTmlz6pQpU2TFihXy7bffyqBBgyQ6OlomT55coe+vysGtDXpxhMov/lDf71DcZ8V+hw7OdejgXJevqr/rQybUffTRR3ZMumeffdZ2jvj+++/lqquukpycnEBvGgAAQLWFTKhTr776qp0AAADcJiSGNHGjw4cPy6hRo+zXUBKK+x2K+6zY79A535xrzrXbHa6l3+NhOg5vjf4EAAAA1Dhq6gAAAFyAUAcAAOAChDoAAAAXINQ5UHJysh1LT587q8+nnTlzppx99tnlfk///v3FGOMzFRQUSDAZOXJkiX3YsGFDud9z88032zK6r2vXrpWrr75ags2WLVtK7LdOZQ2MHYznunv37jJnzhw7+KZub1JSUokyzzzzjB0P8tChQ7Jo0SJp3br1Cd/34YcftsdP9/+bb76RCy64QIJlv/V51KNHj7bX7a+//mrL6FiaCQkJfv+cOOlc69igxbf/008/dfW5VqV9xnV6/PHHg/ZcV+RvVd26de3vMn2yk4699vHHH0vjxo1P+N5V+X3ghH1u0KCBvPLKK7Jx40a77du2bZOXX35ZYmNjy33fqn4uiiPUOVDPnj3t0CsXX3yx9O3bVyIjI2XhwoVyyimnlPt9Bw4csGPweabmzZtLsFm3bp3PPlx66aVllu3atat88MEHMmnSJOnUqZPMmjXLTm3btpVgon+ciu5znz597PJ//OMfrjnXOtD3mjVr5JFHHil1/bBhw+TRRx+VgQMHykUXXSQHDx6UBQsW2D8IZbn11lvtoOL6y79z5872/fV7dDzKYNhv/Tzrdj/33HP264033ijnnHOODQb+/Jw47Vwr/WNVdPvvuOOOct8z2M+1Krq/Og0YMEAKCwtl+vTpQXuuK/K3Ki0tTa677jq55ZZbbPkzzjhDZsyYUe77VuX3gVP2+YwzzrCThvV27drJvffea8fE1b9TJ1LZz0VZtPcrk4OPQcOGDY3q3r17mWX69+9v9u/fH/Btrc40cuRIs3r16gqX//DDD80///lPn2Xp6enm9ddfD/i+VGdKS0szP//8s2vPtUpKSvJZtmvXLjN06FDvfGxsrCkoKDC33XZbme/zzTffmPHjx3vnw8LCzL/+9S8zfPjwgO9jRfe7+HT++efbck2bNvXb58Rp+zx58mQzc+bMSr2PG8+1HoPFixeXWyaYznVpf6v0c3z48GFz0003ecucc845tsxFF11U5vtU5feBU/ZZSpluvvlm89tvv5k6deqUWaYqn4vSJmrqgkC9evXs13379pVb7tRTT5WtW7fK9u3bbY3VeeedJ8HmrLPOss0XmzZtkqlTp0rTpk3LralbvHixzzL9b06XByv9r+/uu++Wt99+2/Xn2qNFixa2ybHoudSmjeXLl5d5LvU4denSxed7tLlC54P5/OtnXWtvcnNz/fY5caJevXrZpittonrttdfktNNOK7OsG8+1Nj9ec801Faq9CaZzXfxvlZ63qKgon3P3448/2ibJss5dVX4fOP3vc7169ew+HD9+3G+fi7IQ6hwuLCxMxo0bJ19++aWsX7++zHL6QbnvvvvsfRwaCsLDw+Xrr7+2Dw8OFvqh9VRVP/TQQ/bD/cUXX9gAUxqtntYPQFE6r8uDVb9+/aR+/fryzjvvuPpcF+U5X5U5lw0bNrT3pLnp/GvT0pgxY+wtBeU997GynxOnmT9/vtxzzz3Su3dvGT58uG3O0mYnvY5D5VzrfbF6jk/UDBlM57q0v1V6fnSwXb1dpKLnriq/D5z89/n000+XESNGyJtvvunXz0V5Al59yVT2MXjttdfMli1bTJMmTSp1nCIiImwT3rPPPhu0x7devXomNzfX3HfffaWu12r922+/3WfZQw89ZLKysgK+7VWd5s+fb+bMmePqc128aapr1652WXx8vE+5v//977aJvbT3SEhIsN9z8cUX+ywfM2aMbaoLhv0ufg5nz55tVq5caWJiYvz6OXHqPnumFi1a2HKXX355SJxrnTZs2GBeeeWVSr+vk891aX+r7rjjDtvsWLzs8uXLzejRo0t9n6r8PnDSPkuRST/Leo3OmzfPfsYr894n+lyUNVFT52Djx4+Xa6+9Vi677DJb/V4Zx44dk9WrVzumx1BV6H93P/30U5n7kJWVJXFxcT7LdF6XB6NmzZrZThJvvfVWSJ1rz/mqzLnUnnS63244/1oL9dFHH9nOLnrjdXm1dFX5nDid9mjdvXt3mdvvpnOttKNDmzZtKv05d/K5LutvlZ4frYH2NFFW5NxV5feBE/8+n3rqqbb2TT/PN9xwg72G/fm5KAuhzqH0gtEL4fLLL7f3TlWWVtkmJiZKZmamBCvtTdaqVasy9yE9Pd1WVRelfxR1eTDS3nA5OTkyd+7ckDrX+stLt73ouYyJibG93so6l0ePHpWVK1f6fI82heh8MJ1/T6DT+6Y00J/ovtmqfE6cTm8b0CaqsrbfLefa4/7775cVK1bYoWzccK7L+1ul5+3IkSM+506H/9B/YMo6d1X5feC0v88xMTG2R6zu+/XXX1+l572e6HNRnoBXYTL5HoNXX33V9m7s0aOHiYuL804nnXSSt8yUKVPMCy+84J0fMWKE6du3r62y7dSpk5k2bZo5dOiQOffcc4Pm+P71r3+1+9y8eXNbBb9w4UKTk5NjexeVts9a5siRI2bIkCG2R5X2FNMm2bZt2wZ8Xyo7aW++rVu3mtTU1BLr3HCuo6OjTYcOHeykBg0aZF97enkOGzbM7Nu3z1x33XWmXbt2thfYpk2bTN26db3voT0FH3nkEe/8rbfeanvE3XPPPaZNmzbmf//3f+17NG7cOCj2W5tjZs2aZbZv327at2/v81mPjIwsc79P9Dlx8j7ruhdffNH2fNTt16alFStWmB9//NFERUW59lwXbY779ddfzYMPPljqewTbua7I3yptotTfbb169TKdO3c2X331lZ2KN0f369fPO1+R3wdO3eeYmBg7CsOaNWtMy5YtfcqEh4eXus8V/VxUcAr8hcHkewzKokNZeMosWbLEdoH2zI8dO9Z+cPT+hczMTPPJJ5+Yjh07BtWx/eCDD8zOnTvtPuzYscPO64eirH32dBXfuHGj/Z6MjAxz9dVXB3w/qjJpSFNnnXVWiXVuONc9e/Ys9Zouul/PPPOM3R/9471o0aISx0LvXdHgXnSZ/gH0HAu9d+XCCy8Mmv3WX95l0e8ra79P9Dlx8j7rHz69bzQ7O9v+A6b79sYbb5QIZ247154yDzzwgDl48KAdoqO09wi2c12Rv1UaxCZMmGD27t1rA+306dNtwCn+PkW/pyK/D5y6zz3LuA6UfuZL2+eKfi4qMoX9+wUAAACCGPfUAQAAuAChDgAAwAUIdQAAAC5AqAMAAHABQh0AAIALEOoAAABcgFAHAADgAoQ6AAAAFyDUAUAxxhhJSkriuAAIKoQ6ACFj8uTJNrDppA/bzsrKsg/eHjBggH1IvEd8fLx8+umnFXpPAiAApyDUAQgpGtY0tJ155ply9dVXy5IlS+Tll1+WTz75ROrUqWPLZGdn29AHAMEm4A/IZeIYcA1wDdTGNaAPV585c2aJ5Zdddpl9wPb999/vfdh2UlKSfR0ZGWnGjx9vdu3aZR8urg+VT05Otuv0wdtF6bwu14euz5o1y2RlZZn8/Hzz7bffmt69e/v8TC2bkpJiJk2aZPLy8sy2bdvsA9+LlmnSpImZNm2a92Ho3333nc+D7K+//nqzcuVKu12bNm0yTz/9tKlTpw6fJz5PXAMSsscg4BvAxDHgGuAaCGio02n16tVm7ty5JULd0KFDbeC69NJLTbNmzcwll1xibr/9druuYcOGtmz//v1NXFycndfl7du3N//93/9t2rZta1q3bm2effZZc+jQIdO0aVOfULdnzx7z0EMPmVatWpnhw4ebY8eOmbPPPtuuj46ONr/88otZunSp/Zla5pZbbjEXX3yxXa/bk5uba+655x7TokUL06dPH7N582Yb7Lie+J3CNSChegwCvgFMHAOuAa6BgIe6Dz74wKxfv75EqHv55ZfN4sWLy3zPomXLmzIyMswjjzziE+reffddnzJas/fggw/a11prd+DAAdOgQYNS32/RokXeGkPPdNddd5mdO3fyeeLzxDUgoXkMIgLd9gsATqAdJbTTQ3HvvPOOLFq0SH788UeZP3++vfdO58sTHR0to0aNkmuuuUYSEhIkIiJCTj75ZGnWrJlPubVr1/rMa8eNxo0b29cdO3aU1atXy/79+0v9GR06dJBLLrlEnnzySe8yvSdQf45OBQUFldp/AMGPUAcAInLuuefKli1bShwLDVYtWrSwnSr69OkjH330kSxevFhuueWWMo/b3/72N+nbt688/vjj8ssvv9iA9fHHH0tUVJRPuaNHj/rMa6gMD/+9/9qJQtmpp54qI0eOlBkzZpRY99tvv3FOgRBEqAMQ8i677DJp3769pKWllXos8vPzbZjTScPZggULpEGDBrYWTXvJenrNemgNmtbwzZo1y1tzp71tK0Nr8f74xz96f05xq1atknPOOUc2bdoU8ucPwO8IdQBCSt26dSUuLs4GMf161VVXSUpKivzzn/+Ud999t0T5wYMHS2Zmpq2xKywstDV0Op+bm2vXb926VXr37i1fffWVHD582C7/+eef5cYbb7TvqbVvzz33nLcGrqI++OAD+fOf/2yDoW6f/sxOnTrJrl275JtvvpFnn33WNgVv377dBk3dNm2SbdeunYwYMcJvxwtA8GCcOgAhRZtR9d41DWN6j5zW0j366KP2CRIajEqrpRs2bJisWLFCvvvuO1vj9l//9V/e+++GDh1qm1p37Nhhg58aMmSIrV37+uuvbbDTmj2tWasMbZq94oorJCcnR+bNmycZGRmSnJwsx48ft+t10ORrr73WltHt0qCnAXTbtm1+OU4Ago8OoV7yzmAAAAAEFWrqAAAAXIBQBwAA4AKEOgAAABcg1AEAALgAoQ4AAMAFCHUAAAAuQKgDAABwAUIdAACACxDqAAAAXIBQBwAA4AKEOgAAABcg1AEAAEjw+//mag7y9XgPOgAAAABJRU5ErkJggg==" }, "metadata": {}, "output_type": "display_data", "jetTransient": { "display_id": null } } ], "execution_count": 6 }, { "metadata": {}, "cell_type": "markdown", "source": [ "As remarked by Hubble, the relationship looks pretty linear. As done before us, we will fit a linear regression model to these data points\n", "\n", "## Linear Regression\n", "\n", "The following code fits a linear regression model. As with Hubble, we will fit a model that does *not* include an intercept term. This is done in `delicatessen` by modifying the specified design matrix in `X`. Given that we have relatively few data points, we will also apply a finite-sample correction to the sandwich variance estimator. " ], "id": "a6b9d0c001526f12" }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-10T15:06:54.692202400Z", "start_time": "2026-04-10T15:06:54.655737200Z" } }, "cell_type": "code", "source": [ "def psi_lr(theta):\n", " return ee_regression(theta=theta, X=d[['x']], y=d['y'], model='linear')" ], "id": "45f9deda06dde4b5", "outputs": [], "execution_count": 21 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-10T15:06:55.908530800Z", "start_time": "2026-04-10T15:06:55.879419500Z" } }, "cell_type": "code", "source": [ "estr1 = MEstimator(psi_lr, init=[0., ], finite_correction='HC1')\n", "estr1.estimate()" ], "id": "3e559167e0c3e2ad", "outputs": [], "execution_count": 22 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-10T15:06:57.348597600Z", "start_time": "2026-04-10T15:06:57.248407100Z" } }, "cell_type": "code", "source": "estr1.print_results()", "id": "dfd005693e9b553", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================================\n", " Estimation Method: M-estimator\n", "--------------------------------------------------------------\n", "No. Observations: 24 | No. Parameters: 1\n", "Solving algorithm: lm | Max Iterations: 5000\n", "Solving tolerance: 1e-09 | Allow P-Inverse: 1\n", "Derivative Method: approx | Deriv Approx: 1e-09\n", "Small N Correction: HC1 | Distribution: t-stat\n", "==============================================================\n", " Theta StdErr Z-score LCL UCL P-value S-value \n", "--------------------------------------------------------------\n", " 76.58 4.93 15.54 66.38 86.78 0.00 43.05 \n", "==============================================================\n" ] } ], "execution_count": 23 }, { "metadata": {}, "cell_type": "markdown", "source": [ "This estimate is the Hubble Constant.\n", "\n", "From the Hubble Constant, we can estimate the age of the universe. To do that, we simply take the inverse of the Hubble Constant (and convert units). However, simply inverting the Hubble Constant does not directly give us the variance for the estimated age of the universe. Here, the ability of estimating equations to automate the delta method comes in handy. We can simply stack another estimating equation that applies our transformation.\n", "\n", "The following code (re-)estimates the Hubble Constant along with Hubble Time." ], "id": "cb655fa11fbeea45" }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-10T15:23:05.461956600Z", "start_time": "2026-04-10T15:23:05.385230200Z" } }, "cell_type": "code", "source": [ "def psi_htime(theta):\n", " h_constant = theta[:1]\n", " h_time = theta[1]\n", "\n", " # Estimating function for Hubble Constant\n", " ee_lr = ee_regression(theta=h_constant, X=d[['x']], y=d['y'], model='linear')\n", "\n", " # Estimating function for Hubble Time\n", " # NOTE: the division by 997.8 is to convert time to billions of years\n", " ee_htime = 1 - h_constant[0]*h_time/997.8 * np.ones(d.shape[0])\n", "\n", " # Returning stacked estimating functions\n", " return np.vstack([ee_lr, ee_htime])" ], "id": "ead3864c8be4521d", "outputs": [], "execution_count": 56 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-10T15:23:07.305336200Z", "start_time": "2026-04-10T15:23:07.236863800Z" } }, "cell_type": "code", "source": [ "estr2 = MEstimator(psi_htime, init=[76., 0], finite_correction='HC1')\n", "estr2.estimate()" ], "id": "d9876058a310b5e9", "outputs": [], "execution_count": 57 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-10T15:23:36.665108500Z", "start_time": "2026-04-10T15:23:36.606408Z" } }, "cell_type": "code", "source": "estr2.print_results(decimals=1)", "id": "be07e33c67e3169d", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================================\n", " Estimation Method: M-estimator\n", "--------------------------------------------------------------\n", "No. Observations: 24 | No. Parameters: 2\n", "Solving algorithm: lm | Max Iterations: 5000\n", "Solving tolerance: 1e-09 | Allow P-Inverse: 1\n", "Derivative Method: approx | Deriv Approx: 1e-09\n", "Small N Correction: HC1 | Distribution: t-stat\n", "==============================================================\n", " Theta StdErr Z-score LCL UCL P-value S-value \n", "--------------------------------------------------------------\n", " 76.6 5.0 15.2 66.1 87.0 0.0 41.3 \n", " 13.0 0.9 15.2 11.3 14.8 0.0 41.3 \n", "==============================================================\n" ] } ], "execution_count": 59 }, { "metadata": {}, "cell_type": "markdown", "source": [ "So, the estimated age of the universe is 13.0 billion years (95% CI: 11.3, 14.8).\n", "\n", "## Robust Linear Regression\n", "\n", "Given that we see some at least one outlier in the previous plot, we might be concerned about our estimates from before. Here, we repeat the process of estimating the Hubble Constant and Time but using robust linear regression to dampen the influence of any outliers. The following code applies this process" ], "id": "668b4ea4a03d831a" }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-10T15:26:30.830545400Z", "start_time": "2026-04-10T15:26:30.806807200Z" } }, "cell_type": "code", "source": [ "def psi_rlr(theta):\n", " h_constant = theta[:1]\n", " h_time = theta[1]\n", "\n", " # Estimating function for Hubble Constant\n", " ee_lr = ee_robust_regression(theta=h_constant, X=d[['x']], y=d['y'], model='linear',\n", " loss='huber', k=10.)\n", "\n", " # Estimating function for Hubble Time\n", " # NOTE: the division by 997.8 is to convert time to billions of years\n", " ee_htime = 1 - h_constant[0]*h_time/997.8 * np.ones(d.shape[0])\n", "\n", " # Returning stacked estimating functions\n", " return np.vstack([ee_lr, ee_htime])" ], "id": "e6242cb0c89b0a8b", "outputs": [], "execution_count": 63 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-10T15:26:32.481795600Z", "start_time": "2026-04-10T15:26:32.447110400Z" } }, "cell_type": "code", "source": [ "estr3 = MEstimator(psi_rlr, init=[76., 13., ], finite_correction='HC1')\n", "estr3.estimate()" ], "id": "536bbb7a28ef0420", "outputs": [], "execution_count": 64 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-10T15:26:34.233837Z", "start_time": "2026-04-10T15:26:34.170108700Z" } }, "cell_type": "code", "source": "estr3.print_results()", "id": "9159160b007335bc", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================================\n", " Estimation Method: M-estimator\n", "--------------------------------------------------------------\n", "No. Observations: 24 | No. Parameters: 2\n", "Solving algorithm: lm | Max Iterations: 5000\n", "Solving tolerance: 1e-09 | Allow P-Inverse: 1\n", "Derivative Method: approx | Deriv Approx: 1e-09\n", "Small N Correction: HC1 | Distribution: t-stat\n", "==============================================================\n", " Theta StdErr Z-score LCL UCL P-value S-value \n", "--------------------------------------------------------------\n", " 74.44 4.57 16.29 64.97 83.92 0.00 43.30 \n", " 13.40 0.82 16.29 11.70 15.11 0.00 43.30 \n", "==============================================================\n" ] } ], "execution_count": 65 }, { "metadata": {}, "cell_type": "markdown", "source": [ "Here, we get a slightly lower estimate for the Hubble Constant (74.4) and an estimated age of the universe of 13.4 billion years (95% CI: 11.7, 15.1). This is fairly close to the modern estimates of the age of the universe (13.787 billion years, Planck Collaboration 2020). It is also well within the random error of our model.\n", "\n", "## References\n", "\n", "Planck Collaboration (2020). Planck 2018 results. VI. Cosmological parameters. *Astron. Astrophys*, 641, A6.\n", "\n", "Hubble E. (1929). A relation between distance and radial velocity among extra-galactic nebulae. *Proceedings of the National Academy of Sciences*, 15(3), 168-173." ], "id": "60adcd121eb5ec8d" } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.6" } }, "nbformat": 4, "nbformat_minor": 5 }