{ "cells": [ { "metadata": {}, "cell_type": "markdown", "source": [ "# Michelson (1880): Speed of Light\n", "\n", "Michelson, as part of the US Naval Academy, sought to measure the (two-way) speed of light using a revolving mirror. The experimental setup can be found in the associated reference. Here, we apply estimating equations to re-analyze this historic data. This data can be found from a variety of online resources, including the R package `morley`.\n", "\n", "## Setup" ], "id": "3fac94f1c79cc1fc" }, { "cell_type": "code", "id": "initial_id", "metadata": { "collapsed": true, "ExecuteTime": { "end_time": "2026-04-29T18:27:57.152289400Z", "start_time": "2026-04-29T18:27:56.022399900Z" } }, "source": [ "import numpy as np\n", "import scipy as sp\n", "import pandas as pd\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt\n", "\n", "import delicatessen\n", "from delicatessen import MEstimator\n", "from delicatessen.estimating_equations import ee_mean_robust\n", "\n", "print(\"Versions\")\n", "print(\"NumPy: \", np.__version__)\n", "print(\"SciPy: \", sp.__version__)\n", "print(\"Pandas: \", pd.__version__)\n", "print(\"Matplotlib: \", mpl.__version__)\n", "print(\"Delicatessen:\", delicatessen.__version__)" ], "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Versions\n", "NumPy: 2.3.5\n", "SciPy: 1.16.3\n", "Pandas: 2.3.3\n", "Matplotlib: 3.10.8\n", "Delicatessen: 4.2\n" ] } ], "execution_count": 1 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-29T18:27:57.171189200Z", "start_time": "2026-04-29T18:27:57.152289400Z" } }, "cell_type": "code", "source": [ "d = pd.read_csv(\"data/speed1879.csv\")\n", "\n", "# The copy of data comes from MASS, which has different units,\n", "# so we convert them to megameters per second\n", "d['Speed'] = (d['Speed'] + 299000) / 1000" ], "id": "747e36318b102cdb", "outputs": [], "execution_count": 2 }, { "metadata": {}, "cell_type": "markdown", "source": "We can begin by plotting the speed of light data. Here, the units are megameters (1000 km) per second", "id": "3cbf1f94f421fbc0" }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-29T18:27:57.325397200Z", "start_time": "2026-04-29T18:27:57.172381700Z" } }, "cell_type": "code", "source": [ "plt.hist(d['Speed'], bins=20)\n", "plt.xlabel(\"Speed (1000 km/s)\")\n", "plt.tight_layout()" ], "id": "b1061fb2ad780f85", "outputs": [ { "data": { "text/plain": [ "
" ], "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAKRBJREFUeJzt3QmUFdWBP+DL1ihMgxoFgbgh4DIoyqKgoBhGguOCS9znuEVHlLgFJ8jonMG1QaMYEaPjApqoR0dFHVEwGDGi4C6ioGhEJAgogrIpi9b/3Jrp/ndDg2C6ed2X7zvnntdVdd+r6rrV3b++7956dUIIWQAAoNarW+gDAACgagh2AACJEOwAABIh2AEAJEKwAwBIhGAHAJAIwQ4AIBGCHQBAIuqHGqhly5ZhyZIlhT4MAIAaobi4OHz22We1L9jFUDdnzpxCHwYAQI3SqlWrHwx3NS7YlfbUxYPXawcAbO6Ki4vzTq8NyUU1LtiVigcv2AEAbDiTJwAAEiHYAQAkQrADAEiEYAcAkAjBDgAgEYIdAEAiBDsAgEQIdgAAiRDsAAASIdgBACRCsAMASIRgBwCQCMEOACARgh0AQCIEOwCARAh2AACJqF/oAwDScePUSZtkPwP26rZJ9gNQ2+ixAwBIhGAHAJAIwQ4AIBGCHQBAIgQ7AIBECHYAAIkQ7AAAEiHYAQAkQrADAEiEYAcAkAjBDgAgEYIdAEAiBDsAgEQIdgAAiRDsAAASIdgBACRCsAMASIRgBwCwuQa7Hj16hCeffDLMmTMnZFkW+vbtu1ad3XffPTzxxBPhq6++CkuXLg2vvvpq2GGHHarqmAEAqIpg17hx4zBlypTQv3//Sre3bt06TJw4Mbz//vuhZ8+eYe+99w5XX311+Pbbbzd2VwAAbIT6YSONHTs2L+ty7bXXhqeffjoMHDiwbN3HH3+8sbsBAKCQY+zq1KkTDj/88DBjxow8/M2fPz9Mnjy50rdrSxUVFYXi4uIKBQCAAge7Zs2a5cHssssuy4Nd7969w+jRo8Njjz0WDjrooEqfM2jQoLB48eKyEsfuAQBQ4GBXt+7/vlycOHHzzTfnY/GGDh0annrqqdCvX79Kn1NSUhKaNGlSVlq1alWVhwQAsNnY6DF267NgwYKwatWqMG3atArrp0+fHrp3717pc1auXJkXAABqUI9dDHWvvfZa2G233Sqsb9euXZg1a1ZV7goAgL+3xy7e7qRNmzZly7vsskvo0KFDWLhwYZg9e3a44YYbwkMPPRT+8pe/hOeffz706dMnHHnkkfmtTwAAqEHBrnPnzmHChAlly8OGDcsfR40aFc4888zw+OOP5+Pp4qSIW265JXzwwQfhuOOOCy+99FLVHjkAABXUCSFkoQaJs2rj7Ng4kWLJkiWFPhxgI9w4ddImOV8D9uq2SfYDUNuykc+KBQBIhGAHAJAIwQ4AIBGCHQBAIgQ7AIBECHYAAIkQ7AAAEiHYAQAkQrADAEiEYAcAkAjBDgAgEYIdAEAiBDsAgEQIdgAAiRDsAAASIdgBACRCsAMASIRgBwCQCMEOACARgh0AQCIEOwCARAh2AACJEOwAABIh2AEAJEKwAwBIhGAHAJAIwQ4AIBGCHQBAIgQ7AIBECHYAAIkQ7AAAEiHYAQAkQrADAEiEYAcAsLkGux49eoQnn3wyzJkzJ2RZFvr27bvOur///e/zOhdddNHfe5wAAFR1sGvcuHGYMmVK6N+//3rrHX300aFr1655AAQAoPrV39gnjB07Ni/r07JlyzB8+PDw85//PIwZM+bvOT4AAAo1xq5OnTrhD3/4Q7jhhhvCtGnTqvrlAQCoqh67HzJw4MCwevXqcMstt2xQ/aKiotCwYcOy5eLi4qo+JACAzUKVBruOHTvmEyXi44YaNGhQGDx4cFUeBkCtc+PUSZtsXwP26rbJ9gXU4rdi44zZZs2ahU8//TSsWrUqLzvvvHO48cYbw8yZMyt9TklJSWjSpElZadWqVVUeEgDAZqNKe+zi2Lrx48dXWDdu3Lh8/ciRIyt9zsqVK/MCAMAmDnbxdidt2rQpW95ll11Chw4dwsKFC8Ps2bPzx/Jir928efPCjBkz/s5DBQCgSoNd586dw4QJE8qWhw0blj+OGjUqnHnmmRv7cgAAFCrYvfDCC/ktTTZU7NEDAKD6+axYAIBECHYAAIkQ7AAAEiHYAQAkQrADAEiEYAcAkAjBDgAgEYIdAEAiBDsAgEQIdgAAiRDsAAASIdgBACRCsAMASIRgBwCQCMEOACARgh0AQCIEOwCARAh2AACJEOwAABIh2AEAJEKwAwBIhGAHAJAIwQ4AIBGCHQBAIgQ7AIBECHYAAIkQ7AAAEiHYAQAkQrADAEiEYAcAkAjBDgAgEYIdAEAiBDsAgEQIdgAAiRDsAAA212DXo0eP8OSTT4Y5c+aELMtC3759y7bVr18/DBkyJLzzzjth6dKleZ177703tGjRoqqPGwCAvzfYNW7cOEyZMiX0799/rW2NGjUKHTt2DFdffXX+eOyxx4bddtstD4IAAFSv+hv7hLFjx+alMosXLw69e/eusO5Xv/pVeO2118IOO+wQZs+e/eOPFACAqg12G6tp06bh+++/D1999VWl24uKikLDhg3LlouLi6v7kAAAklStkydiYBs6dGh48MEHw5IlSyqtM2jQoLynr7TEcXkAANSgYBcnUjz88MOhTp064bzzzltnvZKSktCkSZOy0qpVq+o6JACApNWvzlC30047hZ/97Gfr7K2LVq5cmRcAAGpYsCsNdW3btg2HHHJIWLhwYVXvAgCAqgh28XYnbdq0KVveZZddQocOHfIAN3fu3PDII4/ktzo54ogjQr169ULz5s3zenH7qlWrNnZ3AABUV7Dr3LlzmDBhQtnysGHD8sdRo0aFwYMHl92wON7rrryePXuGF154YWN3BwBAdQW7GM7ihIh1Wd82AACqj8+KBQBIhGAHAJAIwQ4AIBGCHQBAIgQ7AIBECHYAAIkQ7AAAEiHYAQAkQrADAEiEYAcAkAjBDgAgEYIdAEAiBDsAgEQIdgAAiRDsAAASIdgBACSifqEPAKAmu3HqpEIfAsAG02MHAJAIwQ4AIBGCHQBAIgQ7AIBECHYAAIkQ7AAAEiHYAQAkQrADAEiEYAcAkAjBDgAgEYIdAEAiBDsAgEQIdgAAiRDsAAASIdgBACRCsAMASIRgBwCwuQa7Hj16hCeffDLMmTMnZFkW+vbtu1adK6+8Mnz22Wdh+fLl4U9/+lNo06ZNVR0vAABVFewaN24cpkyZEvr371/p9t/85jfhwgsvDP369Qv7779/WLZsWRg3blxo2LDhxu4KAICNUD9spLFjx+ZlXS6++OJwzTXX5L160WmnnRbmz58fjj766PDQQw9t7O4AACjEGLtddtkltGjRIowfP75s3eLFi8Mrr7wSunXrVpW7AgDg7+2xW5/tt98+f4w9dOXF5dJtayoqKqrwNm1xcXFVHhIAwGajSoPdjzFo0KAwePDgQh8GUIvcOHVSoQ8BIP23YufNm5c/Nm/evML6uFy6bU0lJSWhSZMmZaVVq1ZVeUgAAJuNKg12M2fODHPnzg29evWq8NZqnB07aVLl/2GvXLkyLFmypEIBAGATvBUbb3dS/r50ccJEhw4dwsKFC8Ps2bPDzTffHK644orw4Ycf5kHv6quvzu9p9/jjj/+IwwMAoNqCXefOncOECRPKlocNG5Y/jho1Kpx55pnh+uuvz8Pff/3Xf4WtttoqTJw4MfTp0yesWLFiY3cFAMBGqBNCyEINEt+6jbdIiePtvC0LtYtJDbXDgL3cfgpqk43JRj4rFgAgEYIdAEAiBDsAgEQIdgAAiRDsAAASIdgBACRCsAMASIRgBwCQCMEOACARgh0AQCIEOwCARAh2AACJEOwAABIh2AEAJEKwAwBIhGAHAJAIwQ4AIBGCHQBAIgQ7AIBECHYAAIkQ7AAAEiHYAQAkQrADAEiEYAcAkAjBDgAgEYIdAEAiBDsAgEQIdgAAiRDsAAASIdgBACRCsAMASIRgBwCQCMEOACARgh0AQCIEOwCARFR5sKtbt2646qqrwscffxyWL18ePvroo3DFFVdU9W4AAFhD/VDFBg4cGM4777xw+umnh/feey907tw5jBw5Mnz99ddh+PDhVb07AACqK9gdcMAB4YknnghPP/10vjxr1qxw8sknh/3226+qdwUAQHW+Ffvyyy+HXr16hbZt2+bLe++9d+jevXt45plnKq1fVFQUiouLKxQAAGpAj92QIUNCkyZNwvvvvx++++67UK9evXD55ZeHBx54oNL6gwYNCoMHD67qwwAA2OxUeY/dCSecEE499dRwyimnhI4dO+Zj7S699NJw2mmnVVq/pKQkD4KlpVWrVlV9SAAAm4Uq77G74YYb8l67hx56KF9+9913w0477ZT3zN13331r1V+5cmVeAACoYT12jRo1Ct9//32FdfEt2XgbFAAAalGP3f/8z//kY+o+/fTT/HYn++67b/j1r38d7rnnnqreFQAA1RnsLrjggnD11VeH2267LTRr1ix89tln4Y477shvWgwAQC0KdkuXLg2XXHJJXgAA2HQMfAMASIRgBwCQCMEOACARgh0AQCIEOwCARAh2AACJEOwAABIh2AEAJEKwAwBIhGAHAJAIwQ4AIBGCHQBAIgQ7AIBECHYAAIkQ7AAAEiHYAQAkon6hDwCATevGqZM2yX4G7NVtk+wH+P/02AEAJEKwAwBIhGAHAJAIwQ4AIBGCHQBAIgQ7AIBECHYAAIkQ7AAAEiHYAQAkQrADAEiEYAcAkAjBDgAgEYIdAEAiBDsAgEQIdgAAiRDsAAASIdgBACRCsAMASES1BLuWLVuGP/zhD2HBggVh+fLl4Z133gmdOnWqjl0BAPB/6ocqttVWW4WXXnopPP/88+Gwww4LX3zxRWjbtm1YtGhRVe8KAIDqDHYDBw4Ms2fPDmeddVbZuk8++aSqdwMAQHW/FXvUUUeF119/PTz88MNh/vz54c033wxnn332OusXFRWF4uLiCgUAgBrQY9e6detw3nnnhZtuuilcd911oUuXLuGWW24JK1euDPfdd99a9QcNGhQGDx5c1YcBNd6NUydtkv0M2KvbJtkPFFKKP08pfk/Uwh67unXr5r10l19+eXj77bfDnXfemZd+/fpVWr+kpCQ0adKkrLRq1aqqDwkAYLNQ5cFu7ty5Ydq0aRXWTZ8+Pey4446V1o89eUuWLKlQAACoAcEuzojdbbfdKqxr165dmDVrVlXvCgCA6gx2w4YNC127ds3Hzu26667h5JNPDv/6r/8aRowYUdW7AgCgOoNdnBF7zDHH5IHu3XffDf/xH/8RLr744vDAAw9U9a4AAKjOWbHRmDFj8gIAwKbjs2IBABIh2AEAJEKwAwBIhGAHAJAIwQ4AIBGCHQBAIgQ7AIBECHYAAIkQ7AAAEiHYAQAkQrADAEiEYAcAkAjBDgAgEYIdAEAiBDsAgEQIdgAAiahf6AMAqteNUyc5xeDnic2EHjsAgEQIdgAAiRDsAAASIdgBACRCsAMASIRgBwCQCMEOACARgh0AQCIEOwCARAh2AACJEOwAABIh2AEAJEKwAwBIhGAHAJAIwQ4AIBGCHQBAIgQ7AIBEVHuwGzhwYMiyLAwbNqy6dwUAsFmr1mDXuXPncO6554YpU6ZU524AAKjOYNe4ceNw//33h3POOScsWrTIyQYAqK3BbsSIEWHMmDHhueeeq65dAABQTv1QDU488cTQsWPH0KVLlx+sW1RUFBo2bFi2XFxcXB2HBACQvCoPdj/96U/D7373u3DooYeGFStW/GD9QYMGhcGDB1f1YcCPcuPUSc4cALVWlb8V26lTp9C8efPw5ptvhlWrVuWlZ8+e4cILL8y/rlu34i5LSkpCkyZNykqrVq2q+pAAADYLVd5jF8fUtW/fvsK6kSNHhvfffz8MHTo0fP/99xW2rVy5Mi8AANSwYLd06dLw3nvvVVi3bNmy8OWXX661HgCAquOTJwAAElEts2LXdMghh2yK3QAAbNb02AEAJEKwAwBIhGAHAJAIwQ4AIBGCHQBAIgQ7AIBECHYAAIkQ7AAAEiHYAQAkQrADAEiEYAcAkAjBDgAgEYIdAEAiBDsAgEQIdgAAiRDsAAASUb/QBwAb4sapk5woqGX83MKmp8cOACARgh0AQCIEOwCARAh2AACJEOwAABIh2AEAJEKwAwBIhGAHAJAIwQ4AIBGCHQBAIgQ7AIBECHYAAIkQ7AAAEiHYAQAkQrADAEiEYAcAkAjBDgAgEYIdAEAiqjzYXXbZZeHVV18NixcvDvPnzw+jR48O7dq1q+rdAABQ3cHu4IMPDiNGjAhdu3YNhx56aGjQoEF49tlnQ6NGjap6VwAAlFM/VLHDDjuswvIZZ5wRvvjii9CpU6fw4osvVvXuAACormC3pqZNm+aPCxcurHR7UVFRaNiwYdlycXFxdR8SAECSqnXyRJ06dcLNN98cJk6cGN57771K6wwaNCgfj1da5syZU52HBACQrGoNdnGsXfv27cNJJ520zjolJSWhSZMmZaVVq1bVeUgAAMmqtrdihw8fHo444ohw0EEHrbcXbuXKlXkBAKAGBrsY6o455pjQs2fP8Mknn1THLgAAqO5gF99+PeWUU0Lfvn3DkiVLQvPmzfP1X3/9dfj222+rencAAFTXGLvzzz8/bLXVVuGFF14I8+bNKysnnnhiVe8KAIDq7LGLM2EBANj0fFYsAEAiBDsAgEQIdgAAiRDsAAASIdgBACRCsAMASIRgBwCQCMEOACARgh0AQCIEOwCARAh2AACJEOwAABIh2AEAJEKwAwBIhGAHAJAIwQ4AIBH1w2buxqmTQmoG7NWt0IcAAJvF3/YBNexvrh47AIBECHYAAIkQ7AAAEiHYAQAkQrADAEiEYAcAkAjBDgAgEYIdAEAiBDsAgEQIdgAAiRDsAAASIdgBACRCsAMASIRgBwCQCMEOACARgh0AQCIEOwCARFRbsDv//PPDzJkzwzfffBMmT54cunTpUl27AgCguoLdCSecEG666aZw5ZVXho4dO4YpU6aEcePGhe22285JBwCoTcHu17/+dbjzzjvDqFGjwvTp00O/fv3C8uXLw1lnnVUduwMAoDqCXYMGDUKnTp3C+PHjy9ZlWZYvd+vWzUkHAKgm9av6BbfddttQv379MH/+/Arr4/Luu+++Vv2ioqLQsGHDsuXi4uIKj9WtqG69kJpNde42pRTbCaAm8Dej5p+/jdlHlQe7jTVo0KAwePDgtdbPmTOnIMeTgl8tXlzoQwCglvA3o/acvxjwlixZsmmD3YIFC8Lq1atD8+bNK6yPy/PmzVurfklJST7RorxtttkmLFy4MP8GYsBr1arVD34j1B7aNU3aNT3aNE3atfa222efffaD9ao82K1atSq88cYboVevXuGJJ57I19WpUydfvvXWW9eqv3LlyryUt2aIi8uCXXq0a5q0a3q0aZq0a+2yoTmoWt6KjT1w9957b3j99dfDq6++Gi6++OLQuHHjMHLkyOrYHQAA1RXsHn744fyedVdddVXYfvvtw9tvvx369OkTPv/8cycdAKAaZTW1FBUVZf/5n/+ZPxb6WBTt6hrw87q5XQN+Bxe+DbRr4c9XqGWlzv99AQBALVdtnxULAMCmJdgBACRCsAMASES1BbvLLrssv9XJ4sWL848TGz16dGjXrl2FOq1btw6PPfZYPlv266+/Dg899FBo1qxZhTr77rtvePbZZ8OiRYvymx/fcccd+a1Tfkj8+LJ4H72vvvoqLF26ND+WHXbYocq/z81JIds0ft5wZeXSSy+tlu91c1LIdo3bhw8fHmbPnh2WL18e3nvvvXDuuedWy/e5uSlku8bXiLe3ijeYX7ZsWXjmmWdCmzZtquX73Nz069cvTJkyJW+vWF5++eX8rhOl4kd0xnvGxraK9z175JFH1mrT+LfwqaeeytsmXhvXX399qFdv/R/buPXWW4c//vGP+T7jtXDXXXdt0N9iCqNaZmU888wz2emnn57tueee2d5775099dRT2SeffJI1atQo3x4fP/roo+zRRx/N2rdvn5fRo0dnr7zySlanTp28TosWLbIvv/wyu+2227J27dplnTt3ziZOnJj993//93r33bp162zBggXZ0KFDs3322SdfPvLII7Ptttuu4LNVanMpZJs2b968QjnjjDOy7777Lttll10Kfl5qeylku95xxx3Zhx9+mB188MHZTjvtlJ1zzjnZqlWr8p/XQp+X2l4K2a4vv/xy9sILL+T14/Nuv/32CvtWfvw5OOKII7LDDjssa9OmTda2bdvsmmuuyVasWJG3c9we22rWrFnZIYccknXs2DFvi9hmpc+vW7du9s4772TPPvts1qFDh6xPnz7Z559/nl177bXr3e/TTz+dvfXWW9l+++2XHXjggdmMGTOy+++/X1uGGnk9b5odbbvttlnUo0ePfPnQQw/NVq9enRUXF5fVadKkSf7HulevXvly/CU/b968sl8yscRfPtGuu+66zn09+OCD2X333VfoE5t82ZRtumaJf4DGjx9f8HOQYtmU7Tp16tTsiiuuqLDu9ddfz66++uqCn4fUyqZq1xg2otKgEUt8/vz587Nf/vKXBT8PKZYYvs8666y8/WLIO+6448q27bbbbnl77L///vlyDHKx3Zs1a1ZW59xzz82++uqrrEGDBpW+/u67756/RqdOncrW/fznP8+vlRj+C/39K6HCOdhkY+yaNm2aP8bPgC3tLo5vpa1YsaKszrfffhu+//770L1797I68ePGYr1S33zzTf5YWmdN8ePLDj/88DBjxowwduzYvJt58uTJoW/fvtX6/W2ONlWbrim+rRDb+O67767S74dN367xbaSjjjoqtGzZMl/u2bNn/nZhfOuP2tmu8Tmlr1WqdD8b+jPOhqlbt2448cQT87dEJ02aFDp16hSKiorC+PHjy+p88MEHYdasWaFbt275cnycOnVqhQ8MGDduXH59/OM//mOl+4nPiW+/xo8LLRX3Ea+V/fffX3PVMJsk2MWwdfPNN4eJEyfmY2iiGLbi+/tDhw4NW265ZWjUqFH47W9/G+rXrx9atGiR1/nzn/+cf3JFHEfVoEGDsNVWW4UhQ4bk20rrVPZHP35QbhxfEoNd796987ElcRzJQQcdtCm+3c3CpmzTNZ1++un52JHYptTudr3gggvCtGnT8rFYMUDEn9n+/fuHF198UdPW0nZ9//338yBRUlKS14/P+81vfpOP69rQn3HWr3379vnvwBiWb7/99nDMMceE6dOn520V18VxcOXFDo64LYqPcXnN7aXbKhPXr/nJUd99913+T8K6nkPiwW7EiBH5hXjSSSeVrYsDO48//vhw5JFH5pMb4oUYfwnE/wjifwFR/IUf/4gPGDAgH1g9b968MHPmzPyxtM5a31Dd//2W4sSJ+IssDjKNv7jiQNE46JTa16ZrOuuss8L9999foaeB2tmuMdh17do1f+3Y2xCfH4+hV69emrSWtuvq1avDsccem/e8xl6e+LxDDjkkPP300xv8M876xV64ffbZJ+8t+/3vf59/Nvsee+zhtFGmWt+fHj58ePbpp59mO++88zrr/OQnP8maNm2afz137tzs0ksvXatOHA/QuHHjfPBtHB/wi1/8otLXimMEVq5cmV1++eUV1g8ZMqTCAFKl9rRp+dK9e/d8rEccDK4Na/fP6hZbbJGPB/rnf/7nCuvvvPPOfOC/9q2d7Vq+xDFfcWxf/Hry5MnZrbfeql2r4e/Pn/70p3yCSpwwEZW2ZWmJE1cuvvji/Osrr7wynwRRfnu8NqI42bCy1z/zzDOzhQsXVlhXr169fKLT0UcfrU1DjcsU1fsL5W9/+1s+e2dD6seLMg7GjLOo1lUnXmBLly5d68ItX1566aW1Jk889thjZvDU4jYtLSNHjsxee+21Qv/QJFcK0a5x0H4UB3OXXx//QI0bN67g5ySFUuif19IS9x/DYJywUehzkmJ57rnn8t+NpZMnjj322LJtsS0rmzxR/i4RcZJMnDyxrs9lL508EWfZlq6LbWnyRKippXpeeMSIEdmiRYuygw46qMJtKuJ/6aV14i0r4sUWb0dy6qmn5rco+e1vf1vhdfr375/tu++++Uyr888/P1u2bFl2wQUXVKgzffr0Cv81xK/jxX322WfnM7fia8T/LOIU7RpwwmttKWSblgaB+AclzuAq9LlIqRSyXZ9//vl8Zmy83UnsNYi351i+fHnWr1+/gp+X2l4K2a6xNy+2abwd0VFHHZXNnDkze+SRRwp+TlIo1113XT6zOd4eKM5QjssxYP3TP/1T2e1OYg9dz5498yAWOzpiWfN2J2PHjs3f+ejdu3c+Y7n87U66dOmSt2nLli0r3O7kjTfeyLcdcMAB2QcffKCzJNTYUj0vvC7xF3dpnZKSkrzbP4aweJFccskla73Ovffem/+y+fbbb7O33347+5d/+ZdK91X+dUv/q4z32Yl/JGK3c/zlUgNOdq0uhW7T+F9l/KMS/yst9LlIqRSyXWPQuOeee/JepfizGv+YVPbaSu1q1xj84tu/8XVjyLjqqqvWeSsNZePOwV133ZUH5dgeMZDFt2FLQ10sDRs2zN/yjrdAif8Ix/sUxp+z8q+x4447ZmPGjMl/n8Z72N1www35W6ul22Moj2J4LF239dZb50Fu8eLFee/e3Xffnb81r/1CjTsHdf7vCwAAajmfFQsAkAjBDgAgEYIdAEAiBDsAgEQIdgAAiRDsAAASIdgBACRCsAMASIRgB9RKO+20U/zIg9ChQ4f11mvXrl2YO3du+Id/+IdQG8ycOTNcdNFFBT2Gc889Nzz55JMFPQbgxxHsgHXadtttw2233RZmzZoVvv322zwgjR07NhxwwAG15qyVlJSE4cOHh6VLl+bLDRs2DCNHjgzvvPNOWLVqVRg9enSlzzv44IPDG2+8kX/fH374YTj99NPXqnP++efnQeybb74JkydPDl26dKmwPe7r1ltvDQsWLAhLliwJjzzySGjWrFkolB133DEsX748NG7ceL317rnnntCxY8fQvXv3TXZsQNUQ7IB1evTRR8O+++6bh5rY83XUUUeFCRMmhJ/85Ce14qztsMMO4YgjjgijRo0qW1evXr08iN1yyy1h/PjxlT5v5513DmPGjAnPP/982GeffcLNN98c7rrrrtC7d++yOieccEK46aabwpVXXpmHoClTpoRx48aF7bbbrqzOsGHDwpFHHhmOP/74PCi2bNkyPPbYY6FQ+vbtm39Py5YtW2+9GHgfeOCBcOGFF26yYwOqTsE/sFZxDlwDNe8aaNq0af5B4AcddNB660X9+vXLnn766Wz58uXZX//61+y4446rUOenP/1p9tBDD2WLFi3KP5z88ccfr/AB47H88pe/zKZNm5Z988032fTp07PzzjuvwvYuXbpkb775Zr79tddey44++uh83x06dFjnsQ0YMCB79dVX17l95MiR2ejRo9daP2TIkGzq1KkV1j344IPZM888U7Y8efLkbPjw4f//g7fr1Mn+9re/ZQMHDsyXmzRpkq1YsaLCudhtt93yY95///3XeUzxA94vuuiiCuclnref/exn+fLzzz+f3XLLLdmwYcOyhQsXZvPmzcvOPvvsrFGjRtk999yTf0j7hx9+mPXp02et1x4/fnx27rnnln3Q+yuvvJJ/UHx8/YkTJ+YfDl9at0ePHvkHzW+xxRYFvxYV58A1EDbmHLhgnAPXgGtg7WugXr16eUi46aabsqKiovUGuy+++CIPIG3bts2uuuqqbNWqVdnuu++eb69fv3723nvvZXfddVfWvn37fP0f//jHPLw1aNAgr3PKKadkc+bMyY455phs5513zh8XLFiQnXbaafn2xo0bZ/Pnz8+ft+eee2aHH3549tFHH/1gsIsB8rbbbtvoYPfCCy/kwan8ujPOOCP76quv8q/jccfvsW/fvhXqjBo1Kt9n/PqQQw7Jjy8G5PJ1Pvnkk+ziiy/eoGD3b//2b/m5jaG2dHsMdl9//XV2+eWXZ23atMkf47GMGTMmD3hx3YgRI/LnbbnllmXPi8cRg1qLFi3yto1h7vrrr89at26dt0k81zvssENZ/fjc1atX5wHQz4ffka6BUJvOQcEPQHEOXAM19Bo49thj8x622BMXe3SuvfbabK+99qpQJ1ozPE2aNCkPF/HrU089NQ9x5bfHYLRs2bLs0EMPzZdjD9NJJ51UoU4MLC+99FL+9TnnnJMHlYYNG5Ztjz1PPxTs3nrrreyKK67Y6GD3wQcfZJdddlmFdYcddli+v9iDFcNR1LVr1wp1hg4dmvfkxa9PPvnkPEit+dqxlyz2CP5QsIt1YtiNQbb89hjs/vKXv5Qt161bN1uyZEl27733lq1r3rz5Wj2D8XhKey+33nrrDeqNjW1fGq4V58A1EGrFOTDGDlinOB4sjguLY+vipImePXuGN998c62JBJMmTVpreY899si/jrNW27Rpk08eKC0LFy4MW2yxRdh1111Do0aN8u133313hTpXXHFFvj2KrxUnO6xYsWKd+6zMlltumU9+qG0GDBgQzjnnnHzywrRp09baHs9Fqe+//z58+eWXYerUqWXr5s+fnz+Wn6gRx9eVznRdtGhRPoEkjgmM6+JYuu23336t/cSxiLF9gNpDsAPWK4apOMngmmuuCQceeGA+ESFOGNhQ8TYjcXZpnIRQvsTJGHGAfultSGKQKb+9ffv2oWvXrn9X68TZqFtvvfVGP2/evHmhefPmFdbF5a+//joPivF1V69eXWmd+NzS14izYps2bbrOOuvy4osv5pM84gSNdU1uKC/e9mXNdVHduv/7K75BgwahT58+FW5hctZZZ4Vu3bqFl19+OZx44olhxowZYf/996/w/G222SZ88cUX6z1WoGYR7ICNEnuQ1rxdxpoBLC5Pnz49/zr28LVt2zZ8/vnn4a9//WuFsnjx4nz9nDlzQuvWrdfa/sknn+SvEV9r7733zoPSuvZZmbfeeivsueeeG93CsTewV69eFdYdeuihZb2EMUTFsFq+Tp06dfLl0jpx+8qVKyvUiWE23n/vh3obX3311XDYYYeFf//3f8977/5esac19tKV7+mL3n777TBkyJA8sL/77rvhlFNOKdsW2yP2eMZzCNQuBX8/WHEOXAM17xrYZpttsueeey4fIxfH1cVJDb/4xS+yuXPn5hMhSutFn3/+eXbmmWfmkycGDx6cD7rfY489ygbhxzFrf/7zn7Pu3bvnrxMH5P/ud7/LWrVqldeJEy/imLsLLrggf404ySJOVrjkkkvKJk/Efdx3333568bxbjNmzPjBMXZHHHFEPms0jkMrvz6+RnzeE088kR9X/Lr868RjjLNF45i5OJM1ztCNExR69+5dVueEE07IZ+jGMWhx8sHtt9+ez1Jt1qxZWZ049jBOlujZs2fWsWPHfMxg6bjBDZk8ceCBB+YTWMrPko1j7Nac2LHmTNrSdimd3BFn78bzXf77u+666/IxgnEmbBzrGMcwxtnNpXVOP/30fIJKoa9DxTlwDYSNPQcuGufANeAaWPsaiDNh4x//119/PZ9BGYNOnAQRZ72WvwVGFIPPuHHj8qDz8ccfZ8cff3yF14qD+eOM0RjOYp0YGO64446suLi4wuD+eDuTOOEgDtqfMGFCfkuT0u1xIkCcDBG3x3px5uwPBbs4+zPegqR8ICsNQpUpXyeGz9Ljiccbg86ar9+/f/88uMU6cdLEfvvtV2F7nOxx66235t9PPH+PPvpofi42NNiV3nYkTo741a9+9aOD3axZs7JevXqVbYvh87HHHssnZ8Rjj8+PgTzesqW0ztixY8tu3aI4B66BUJvOQcEPQHEOXAO1+BooHyBqYjn//PPzkFLo4yhU2XffffNgHm87s6HPiTNxY09nvBdfoY9fcQ5cA2GjzkH9Qr8PDFCd7rjjjrDVVlvlkzRKP1Zsc1K/fv1wwQUX5JM9NlSLFi3Caaedlo+BBGqXOv+X8AB+lPgO5tFHHx2eeOIJZxCgwAQ7AIBEuN0JAEAiBDsAgEQIdgAAiRDsAAASIdgBACRCsAMASIRgBwCQCMEOACARgh0AQEjD/wOzdrLRo49pogAAAABJRU5ErkJggg==" }, "metadata": {}, "output_type": "display_data", "jetTransient": { "display_id": null } } ], "execution_count": 3 }, { "metadata": {}, "cell_type": "markdown", "source": "So we some variable measurements. From these observations, we will estimate the central location. This can be done using the mean of the measurements, which can be manually computed using `delicatessen` as follows", "id": "cdc666a92ec0318c" }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-29T18:27:57.346168400Z", "start_time": "2026-04-29T18:27:57.329997600Z" } }, "cell_type": "code", "source": [ "def psi_mean(theta):\n", " return d['Speed'] - theta[0]" ], "id": "8eea651dfbe6e089", "outputs": [], "execution_count": 4 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-29T18:27:57.391881700Z", "start_time": "2026-04-29T18:27:57.350070Z" } }, "cell_type": "code", "source": [ "estr = MEstimator(psi_mean, init=[300, ])\n", "estr.estimate()\n", "estr.print_results()" ], "id": "63c2b5b6e3955a7", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================================\n", " Estimation Method: M-estimator\n", "--------------------------------------------------------------\n", "No. Observations: 100 | 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: None | Distribution: Z-stat\n", "==============================================================\n", " Theta StdErr Z-score LCL UCL P-value S-value \n", "--------------------------------------------------------------\n", " 299.85 0.01 38142.12 299.84 299.87 0.00 inf \n", "==============================================================\n" ] } ], "execution_count": 5 }, { "metadata": {}, "cell_type": "markdown", "source": [ "Therefore, we estimate the speed of light to be 299,850,000 meters per second (95% confidence interval: 299840000, 299870000). The current agreed upon speed of light is 299,792,458 m/s, so the estimate here isn't too far off.\n", "\n", "A concern we might have from the previous plot is that there may be some undue dependence on outliers. Rather than try to remove outliers from the data, we can dampen their influence instead. This can be done by using a version of the mean robust to outliers. Here, we use Huber's proposed method. This can be done using the `ee_mean_robust` function in `delicatessen`" ], "id": "1785c688d69d9672" }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-29T18:27:57.426794Z", "start_time": "2026-04-29T18:27:57.391881700Z" } }, "cell_type": "code", "source": [ "def psi_robust(theta):\n", " return ee_mean_robust(theta=theta, y=d['Speed'],\n", " loss='huber', k=0.1)" ], "id": "386bdf464388c1d0", "outputs": [], "execution_count": 6 }, { "metadata": { "ExecuteTime": { "end_time": "2026-04-29T18:27:57.456231300Z", "start_time": "2026-04-29T18:27:57.428332600Z" } }, "cell_type": "code", "source": [ "estr = MEstimator(psi_robust, init=[299.85, ])\n", "estr.estimate()\n", "estr.print_results()" ], "id": "289cef8d519f0415", "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================================\n", " Estimation Method: M-estimator\n", "--------------------------------------------------------------\n", "No. Observations: 100 | 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: None | Distribution: Z-stat\n", "==============================================================\n", " Theta StdErr Z-score LCL UCL P-value S-value \n", "--------------------------------------------------------------\n", " 299.85 0.01 36984.78 299.84 299.87 0.00 inf \n", "==============================================================\n" ] } ], "execution_count": 7 }, { "metadata": {}, "cell_type": "markdown", "source": [ "Here, we see little change when dampening the influence of potential outliers. This suggests that the deviation of the estimated speed of light is not a result of the potential outliers in the data.\n", "\n", "## References\n", "\n", "Michelson, A. A. (1880). Experimental determination of the velocity of light (Vol. 1). US Nautical Almanac Office." ], "id": "76ae72956cc316f" } ], "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 }