{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Analyzing speed of light measurements" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\"\"\n", "\n", "This notebook is an element of the [risk-engineering.org courseware](https://risk-engineering.org/). It can be distributed under the terms of the [Creative Commons Attribution-ShareAlike licence](https://creativecommons.org/licenses/by-sa/4.0/).\n", "\n", "Author: Eric Marsden . \n", "\n", "---\n", "\n", "This notebook contains an introduction to use of Python and SciPy to analyze historical measurements of the speed of light, an important quest in physics. It includes an introduction to the use of **bootstrap resampling** to generate confidence intervals for a measure. See the [associated course materials](https://risk-engineering.org/statistical-modelling/) for background information and to download this content as a Jupyter notebook." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "import numpy\n", "import scipy.stats\n", "import pandas\n", "import matplotlib.pyplot as plt\n", "plt.style.use(\"bmh\")\n", "%config InlineBackend.figure_formats=[\"svg\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our discussion in this notebook is partially based on information from the article by R. J. MacKay and R. W. Oldford, *Scientific Method, Statistical Method and the Speed of Light*, Statistical Science, 2000, 15:3, pp. 254–278 ([available from Project Euclid](https://projecteuclid.org/download/pdf_1/euclid.ss/1009212817)), and on a discussion of successive measurements of the speed of light by S. Prokhovnik and W. Morris, *A review of speed of light measurements since 1676*, CEN Tech. J. 1993, 7:2, pp. 181–183. \n", "\n", "Around 1879, an instructor in physics at the U.S. Naval Academy in Annapolis, A. Michelson, undertook to measure the speed of light, using an adaptation of Foucault's rotating mirror approach. He made around 100 measurements, which you can [download in CSV format](https://risk-engineering.org/static/data/michelson-speed-light.csv)." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Unnamed: 0datedistinctness of imagetemperature (F)position of deflected imageposition of slitdisplacement of image in divisionsdifference between greatest and leastBCorrevolutions per secondradius (ft)value of one turn of screwspeedremarks
00June 5376114.850.300114.550.171.423-0.132257.3628.6720.99614299850Electric light.
11June 7272114.640.074114.560.101.533-0.084257.5228.6550.99614299740P.M. Frame inclined at various angles
22June 7272114.580.074114.500.081.533-0.084257.5228.6470.99614299900P.M. Frame inclined at various angles
33June 727285.910.07485.840.121.533-0.084193.1428.6470.99598300070P.M. Frame inclined at various angles
44June 727285.970.07485.89O.071.533-0.084193.1428.6500.99598299930P.M. Frame inclined at various angles
\n", "
" ], "text/plain": [ " Unnamed: 0 date distinctness of image temperature (F) \\\n", "0 0 June 5 3 76 \n", "1 1 June 7 2 72 \n", "2 2 June 7 2 72 \n", "3 3 June 7 2 72 \n", "4 4 June 7 2 72 \n", "\n", " position of deflected image position of slit \\\n", "0 114.85 0.300 \n", "1 114.64 0.074 \n", "2 114.58 0.074 \n", "3 85.91 0.074 \n", "4 85.97 0.074 \n", "\n", " displacement of image in divisions difference between greatest and least \\\n", "0 114.55 0.17 \n", "1 114.56 0.10 \n", "2 114.50 0.08 \n", "3 85.84 0.12 \n", "4 85.89 O.07 \n", "\n", " B Cor revolutions per second radius (ft) \\\n", "0 1.423 -0.132 257.36 28.672 \n", "1 1.533 -0.084 257.52 28.655 \n", "2 1.533 -0.084 257.52 28.647 \n", "3 1.533 -0.084 193.14 28.647 \n", "4 1.533 -0.084 193.14 28.650 \n", "\n", " value of one turn of screw speed remarks \n", "0 0.99614 299850 Electric light. \n", "1 0.99614 299740 P.M. Frame inclined at various angles \n", "2 0.99614 299900 P.M. Frame inclined at various angles \n", "3 0.99598 300070 P.M. Frame inclined at various angles \n", "4 0.99598 299930 P.M. Frame inclined at various angles " ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mich = pandas.read_csv(\"https://risk-engineering.org/static/data/michelson-speed-light.csv\")\n", "mich.rename(columns={\"velocity of light in air (km/s)\": \"speed\"}, inplace=True)\n", "mich.head()" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2023-07-06T14:33:16.828341\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.6.3, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.xlabel(\"Measured speed (km/s)\")\n", "plt.title(\"Michelson speed of light measurements\")\n", "plt.hist(mich.speed, alpha=0.5);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The mean value from Michelson’s measurements is" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "299852.4" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mich.speed.mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since there is clearly some measurement uncertainty, we would like to report a **95% confidence interval** for the measurements. We can do this using bootstrap resampling, as implemented in the function below." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "def bootstrap_confidence_intervals(data, estimator, percentiles, runs=1000):\n", " replicates = numpy.empty(runs)\n", " for i in range(runs):\n", " replicates[i] = estimator(numpy.random.choice(data, len(data), replace=True))\n", " est = numpy.mean(replicates)\n", " ci = numpy.percentile(replicates, percentiles)\n", " return (est, ci)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(299852.6815, array([299837.295, 299867.705]))" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "est, ci = bootstrap_confidence_intervals(mich.speed, estimator=numpy.mean, percentiles=[2.5, 97.5])\n", "est, ci" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(15.386500000022352, 15.023500000010245)" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "est - ci[0], ci[1] - est" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first return value is the bootstrapped estimate of the mean, and the second return value is the 95% confidence interval. Given that the confidence interval is roughly symmetrical, we could report this measurement as\n", "\n", " 299853 ± 16 km/s (95% CI) \n", "\n", "or as \n", "\n", " 2998523 [299837, 299869] km/s (95% CI). \n", "\n", "Note that a “naked” ± notation for the margin of error is ambiguous, as is a simple interval, because some technical communities use one or two standard deviations instead of a confidence interval. \n", "\n", "At the time, Michelson made a few corrections to his measurements to account for the influence of temperature, and accounting for a number of sources of uncertainty in his measurement method, reported his estimate of the speed of light in a vacuum as\n", "\n", " 299944 ± 51 km/s" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Just three years later in 1882, Simon Newcomb ran another experiment to measure the speed of light. He measured the time required for light to travel from his laboratory on the\n", "Potomac river to a mirror at the base of the Washington Monument and back, a total distance of about 7400 meters. The data are recorded as deviations from 24800 nanoseconds. They were measured on three different days. The “raw” measurements are given below." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "newcomb1 = [28, -44, 29, 30, 26, 27, 22, 23, 33, 16, 24, 29, 24, 40, 21, 31, 34, -2, 25, 19]\n", "newcomb2 = [24, 28, 37, 32, 20, 25, 25, 36, 36, 21, 28, 26, 32, 28, 26, 30, 36, 29, 30, 22]\n", "newcomb3 = [36, 27, 26, 28, 29, 23, 31, 32, 24, 27, 27, 27, 32, 25, 28, 27, 26, 24, 32, 29, 28, 33, 39, 25, 16, 23]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To transform from these units into the speed of light expressed in km/s, we use the following transformation:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "def newcomb_transform(x):\n", " return 1_000_000 * 7.44373 / (x*0.001+24.8)\n", "\n", "newcomb1t = newcomb_transform(numpy.array(newcomb1))\n", "newcomb2t = newcomb_transform(numpy.array(newcomb2))\n", "newcomb3t = newcomb_transform(numpy.array(newcomb3))" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2023-07-06T14:33:17.061029\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.6.3, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(9, 4))\n", "plt.hist(newcomb1t, bins=10, alpha=0.3, label=\"Day 1\")\n", "plt.hist(newcomb2t, bins=10, alpha=0.3, label=\"Day 2\")\n", "plt.hist(newcomb3t, bins=10, alpha=0.3, label=\"Day 3\")\n", "plt.legend()\n", "plt.xlabel(\"Measured speed (km/s)\")\n", "plt.title(\"Newcomb’s speed of light measurements\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the observations from the three days do seem compatible, even if the data from the first day have a larger degree of variance (uncertainty). In fact, there are two **apparent outliers** in the untransformed data (the two negative values). Newcomb decided to discard the measurement of -44, treating it as an outlier (note that there is no magical method for treating outliers in statistical analysis, and the decision to exclude an observation or measurement should not be taken lightly since it may have a significant effect on statistical analyses). \n", "\n", "We also exclude the measurement of -44 and calculate a bootstrap mean and 95% confidence intervals for these observations." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2023-07-06T14:33:17.226995\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.6.3, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "newcomb1.remove(-44)\n", "newcomb_all = numpy.concatenate([newcomb_transform(numpy.array(newcomb1)), \n", " newcomb2t, \n", " newcomb3t])\n", "# a histogram of all non-outlier measurements\n", "plt.hist(newcomb_all, bins=15, alpha=0.5)\n", "plt.title(\"Newcomb’s speed of light measurements (1 outlier excluded)\")\n", "plt.xlabel(\"Measured speed (km/s)\");" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(299820.11243048194, array([299804.11313327, 299839.43270316]))" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "est, ci = bootstrap_confidence_intervals(newcomb_all, numpy.mean, [2.5, 97.5])\n", "est, ci" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(15.999297212983947, 19.320272676297463)" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "est - ci[0], ci[1] - est" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The uncertainty is not quite symmetrical, so we would prefer to report this result as\n", "\n", " 299821 [299804, 299840] km/s (95% CI)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is interesting to note that the defined true value, 299792.458 km/s, is well outside both the confidence interval for both these sets of measurements. (We don’t know of any good explanation for this.)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For more information on these experiments, see the article \n", "\n", "MacKay, R. J.; Oldford, R. W. Scientific Method, Statistical Method and the Speed of Light. Statist. Sci. 15 (2000), no. 3, 254--278. doi:10.1214/ss/1009212817. https://projecteuclid.org/euclid.ss/1009212817" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.11.4" } }, "nbformat": 4, "nbformat_minor": 1 }