![](\"Images/UCERF.png\")
\n",
"\n",
"\n",
"For this activity, you will use the historical seismicity to determine the probability of various-sized earthquakes occurring in specific areas over the next year, and over the next 30 years. The data came from [https://earthquake.usgs.gov/earthquakes/search/](https://earthquake.usgs.gov/earthquakes/search/)\n",
"\n",
"| | San Francisco area (done as class) | Los Angeles area (done individually for your lab report) |\n",
"|--|--|--|\n",
"|Latitude range |36.25 - 38.75 $^{\\circ}$N |33.5-35.5 $^{\\circ}$N |\n",
"|Longitude range |120.75 - 123.25 $^{\\circ}$W|116.75 - 119.75 $^{\\circ}$W|\n",
"|Date range |01/01/1983 - 12/31/2012 |01/01/1983 - 12/31/2012 |\n",
"|Magnitude ranges|2.0-2.9, 3.0-3.9, up to 9.0-9.9|1-1.9, 2.0-2.9, 3.0-3.9, up to 9.0-9.9|\n",
"|Depth range |All |All |\n",
"|Data source |United States Geologic Survey|United States Geologic Survey|\n",
"|Database |[http://neic.usgs.gov](http://neic.usgs.gov) |[http://neic.usgs.gov](http://neic.usgs.gov) |\n"
]
},
{
"cell_type": "markdown",
"id": "4e0d6320-b1b7-4d8a-b570-7de94921c2d7",
"metadata": {},
"source": [
"![](\"Images/Rectangles.png\")
"
]
},
{
"cell_type": "markdown",
"id": "83d789e6-dd66-4af3-b0ed-31a9eeee8969",
"metadata": {
"tags": []
},
"source": [
" Example : San Francisco Area\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "8e60a6ba-30d9-4e65-8bc6-cf68e1d0d947",
"metadata": {},
"source": [
"| Magnitude range | total \\# of earthquakes 1983-2012 (30 years) | average # of earthquakes per year | MRI (mean recurrence interval in years) | One year probability of earthquake occurring | One year probability of earthquake **not** occurring |\n",
"|--|--|--|--|--|--|\n",
"| 2.0-2.9 | 1716 | 57.20 | 0.017 | 1.000 | 0.000 |\n",
"| 3.0-3.9 | 1326 | 44.20 | 0.023 | 1.000 | 0.000 |\n",
"| 4.0-4.9 | 161 | 5.37 | 0.186 | 1.000 | 0.000 |\n",
"| 5.0-5.9 | 13 | 0.43 | 2.308 | 0.433 | 0.567 |\n",
"| 6.0-6.9 | 3 | 0.10 | 10.000 | 0.100 | 0.900 |\n",
"| 7.0-7.9 | 0 | 0.00 | 50 (extrapolated) | 0.020 (extrapolated) | 0.9800 (extrapolated) |\n",
"| 8.0-8.9 | 0 | 0.00 | 300 (extrapolated) | 0.0033 (extrapolated) | 0.9967 (extrapolated) |\n",
"| 9.0-9.9 | 0 | 0.00 | 1600 (extrapolated) | 0.0006 (extrapolated) | 0.9994 (extrapolated) |"
]
},
{
"cell_type": "code",
"execution_count": 37,
"id": "69148f53-5252-43e8-827a-27cc4ced6fe2",
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd"
]
},
{
"cell_type": "code",
"execution_count": 49,
"id": "dd8e9230-b2bf-4ed2-9d6f-f94f7f952db4",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n", "
" ] }, { "cell_type": "markdown", "id": "01abb120-b97b-494a-be2a-3d793561514b", "metadata": {}, "source": [ "**Question 2.2** Note the slight deviation at small magnitudes (~2-3) to fewer than expected detections of small earthquakes. This *kink* is indicative of something called the *magnitude of completeness* of a seismicity catalog, a threshold magnitude below which not all earthquakes may be detected. This bias can be accounted for by excluding earthquakes below the *magnitude of completeness* while estimate the $b$ value (which we will leave to the experts for now!). Why are smaller earthquakes harder to detect? *Hint:* Think signal-to-noise ratio (SNR).\n", "\n", "**Answer:**" ] }, { "cell_type": "markdown", "id": "26fe3d34-5a8f-497e-b884-30b49f42b9fd", "metadata": {}, "source": [ "
\n", "
" ] }, { "cell_type": "markdown", "id": "94769c19-a5a9-441f-9b40-dbaf47f4b143", "metadata": {}, "source": [ "**Question 2.3** Fill the table of joint probabilities of earthquakes happening within a 30-year period for the Los Angeles Area based on the instructions below:" ] }, { "cell_type": "markdown", "id": "aad3be86-93df-4079-a2b6-a8534de8a6a2", "metadata": {}, "source": [ "- First two columns are provided. The values were calculated for you by using the `catalog_filter` function.\n", "- Calculate the average number of earthquakes per year that occurred in each magnitude range.\n", "- Calculate the mean recurrence interval (MRI) for each magnitude range up through magnitude 6.0-6.9. The MRI is defined as the average time between earthquakes. Extrapolate MRI for magnitude 7.0 or above (see hint below).\n", "- Determine the probability of earthquakes of each magnitude range occurring in one year, expressed as either a fractional probability between 0 and 1.0.\n", " - For earthquakes with MRI’s greater than one year: Fractional probability = 1/ MRI and then multiply by 100 to get % probability. (Note that this is equal to the average # of earthquakes per year. But using the 1/ MRI method allows calculation of probabilities for earthquakes that haven’t occurred over the study period - because we have extrapolated MRI’s.)\n", "- Determine the probability of each earthquake magnitude NOT occurring in a year. This will simply be is 1.0 minus the probability of that thing occurring in a year." ] }, { "cell_type": "markdown", "id": "73c7e8fc-8326-4be0-8206-42a1804b40d0", "metadata": {}, "source": [ " \n", "*Hint:* How to extrapolate MRI:\n", "\n", "1. Plot MRI as a function of magnitude to see the relationship. You will find that MRI is an exponential function of magnitude\n", "\n", "2. Fit a straight line $\\log_{10}(\\text{MRI}) = intercept + slope * M$ where $M$ is the magnitude\n", "\n", "3. Use $intercept$ and $slope$ to calculate MRI for magnitude range 7.0-7.9, 8.0-8.9, and 9.0-9.9." ] }, { "cell_type": "code", "execution_count": 4, "id": "cf1582ab-7c78-4639-88b8-10d36312ca81", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/tmp/ipykernel_1728417/4143073906.py:8: RuntimeWarning: divide by zero encountered in true_divide\n", " MRI = 30/n\n" ] }, { "data": { "text/plain": [ "array([3.92223515e-04, 1.39723348e-03, 1.63934426e-02, 1.43540670e-01,\n", " 1.15384615e+00, 1.50000000e+01, inf, inf,\n", " inf])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# the centers of the magnitude ranges\n", "mag = np.array([1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5])\n", "\n", "# number of earthquakes over 30 years\n", "n = np.array([76487,21471,1830,209,26,2,0,0,0])\n", "\n", "# Mean recurrence intervals in years\n", "MRI = 30/n\n", "\n", "# notice that the last 3 elements are inf\n", "MRI" ] }, { "cell_type": "code", "execution_count": 5, "id": "a6492eab-e82d-4adb-9dff-10d40b851584", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAihElEQVR4nO3deXxV5Z3H8c9DNsIawYisBgUiAgISQNC2oWLRupSio0CtuygOVrowFWtHO9MOtky1Oi4VFcENTBEpKm3UTqMoKBDCKgRFIyQBAmhCAtnvb/4gOIgEgrnJOffc7/v18iX3cHLu7yHwfPMs9xxnZoiISPRp4XUBIiLiDQWAiEiUUgCIiEQpBYCISJRSAIiIRKlYrws4lqSkJOvVq5fXZTSZ/fv307p1a6/LaBJBbhuofZEu6O3Lzs7eY2bJxzvP1wHQqVMnVq1a5XUZTSYrK4v09HSvy2gSQW4bqH2RLujtc8591pDzNAUkIhKlFAAiIlFKASAiEqV8GQDOucucc7PKysq8LkVEJLB8GQBm9qqZTWrTpo3XpYiIBJavdwGJiESyRTkFzMzMpbC4nC5JiUwbk8rYwV29LutLCgARkSawKKeA6QvXU15dC0BBcTnTF64H8E0I+HIKSEQk0s3MzP2y8z+kvLqWmZm5HlX0dQoAEZEmUFhcfkLHvaAAEBEJs6J9FcTGuKP+XpekxGaupn4KABGRMNq8cx9jH30PgPiYr3axiXExTBuT6kVZR6UAEBEJk3e27ObKx5dTa8Yrt5/HH648m65JiTiga1IiM8YN8M0CMGgXkIhIWMxbsY17Fm2gT6e2zL4+jc7tE+nftb2vOvwjNVsAOOdOB34FtDezK5vrfUVEmlIoZPwhM5c/v72V9NRkHpl4Dm0SIuNn60ZNATnnZjvnipxzG444fpFzLtc597Fz7i4AM/vEzG5qzPuJiPhJRXUtU+at5s9vb+Wac3vw1LVpEdP5Q+PXAOYAFx1+wDkXAzwKXAycBUxwzp3VyPcREfGVPWWVTHjyff62YSf3XNKX//xBf2JjImtZ1ZlZ4y7gXArwmpn1r3s9ArjPzMbUvZ4OYGYz6l4vONYUkHNuEjAJIDk5eUhGRkaj6vOzsrIygnq/oyC3DdS+SNfY9hWWhXgwu4KSSuPWgQkM6eSvn/pHjRqVbWZpxzuvKaruCmw/7HU+MNw51xH4HTDYOTf9UCAcycxmAbMAUlNTLchP7QnyU4mC3DZQ+yJdY9q3bOse7n8um/jYODImD2VQ96Sw1tacmiIAjvbpBzOzvcBtDbqAc5cBl3Xp0iWshYmINMaC7HymL1xHSsfWzL5+KN07tPK6pEZpigmrfKD7Ya+7AYUncgHdDlpE/MTMeODNLfziL2sZ1rMDCyaPjPjOH5pmBLAS6O2c6wkUAOOBiSdyAY0ARMQvKmtq+eWCdSxaU8hVad343Q8HEBdhi731aew20HnAciDVOZfvnLvJzGqAKUAmsAnIMLONJ3JdjQBExA++2F/Fj59awaI1hUwbk8rvrzg7MJ0/NHIEYGYT6jm+BFjSmGuLiHgpb89+bpizkoLich6eMJjLBwZvRsJfe5fqaApIRLy0Ku9zbnl2FQAv3jyctJQOHlfUNHw5ltEUkIh4ZfHaQiY++QEntYrnldvPC2znDxoBiIgAB3f6PJa1lZmZuQxL6cATPx7CSa3jvS6rSWkEICJRr7o2xC9fXsfMzFzGDurCczcPC3znDz4dAYiINJeS8mpufyGb9z7ey50X9Gbq6N44d/SneQWNAkBEotb2zw9ww5yVfLZ3P3/8l4FcMaSb1yU1K18GgNYARKSprdlezM1zV1JVE+LZG4cz4oyOXpfU7LQGICJRZ9XOGsbPWk5ifAwLbz8vKjt/8OkIQESkKZgZTy39lEfXVDKoRxJPXZtGxzYJXpflGQWAiESFmtoQ9y7eyAsfbGPoqTE8d8u5tIyL8bosT/kyALQGICLhVFZZw5QXV5OVu5vJ6WcwNGFH1Hf+oDUAEQm4HSXlXPn4MpZ+tIcZ4wbwy4vOpEWUbPM8Hl+OAEREwmFDQQk3zV3Jgcpa5twwlG/1Tva6JF9RAIhIIP1j0y7umJdDUmIcCyaPJPXUtl6X5DsKABEJnLnL8vjNqxvp16U9T1+XxintWnpdki/5MgC0CCwi30RtyPjt6x/yzHt5XHhWJx4aP4hW8b7s5nxBi8AiEggHqmq49blsnnkvjxvP68mfrxmizv849KcjIhGvaF8FN81dxcbCEn5zeT+uG5nidUkRQQEgIhFt88593PjMSorLq3nqujS+e2Ynr0uKGAoAEYlY72zZze0vrKZ1QgwZt46gf9f2XpcUURQAIhKR5q3Yxj2LNtD7lDY8c8NQOrdP9LqkiKMAEJGIEgoZf8jM5c9vbyU9NZlHJp5DmwR1Zd+E/tREJGJUVNfys4w1LFm/k2vO7cF9l/UjNsaXmxkjgi8DQJ8DEJEj7Smr5JZnV7FmezH3XNKXm87vGTWPbmwqvoxOfQ5ARA73cVEZP3zsPTbt2MfjPxrCzd86XZ1/GPhyBCAicsjyrXu59blVxMe2YP6kEQzqnuR1SYGhABAR31qQnc/0hetI6dia2dcPpXuHVl6XFCgKABHxHTPjwbc+4uF/fMR5vTry2I+G0D4xzuuyAkcBICK+UllTyy8XrGPRmkKuSuvGb8cOID7Wl8uVEU8BICK+8cX+Km59LpsVeZ8zbUwqt6efocXeJqQAEBHPLMopYGZmLoXF5ZzSLoHakLGvooaHJwzm8oHaBt7UFAAi4olFOQVMX7ie8upaAHbtqwTgJxf0UuffTDSxJiKemJmZ+2Xnf7iXsws8qCY6NdsIwDnXGngMqAKyzOyF5npvEfGfwuLyEzou4deoEYBzbrZzrsg5t+GI4xc553Kdcx875+6qOzwOWGBmtwCXN+Z9RSSyffDJXlrUs7jbJUl39WwujZ0CmgNcdPgB51wM8ChwMXAWMME5dxbQDdhed9rXx30iEniVNbXM+Nsmxj/5PkmtYkk4YntnYlwM08akelRd9HFm1rgLOJcCvGZm/etejwDuM7Mxda+n152aD3xhZq855+ab2fh6rjcJmASQnJw8JCMjo1H1+VlZWRlBvd9RkNsGat83UVAa4s/rKtleGuI73WKZcGY8q4tqeHlLNXsrjI4tHVf0iWNkl6b/wFfQv3+jRo3KNrO0453XFGsAXfn/n/ThYMc/HHgYeMQ5dwnwan1fbGazgFkAqamplp6e3gQl+kNWVhZBbV+Q2wZq34kIhYw5y/K4/4PNtE2I5clrz+bCsw4+tvEi4O6wvMuJCfr3r6GaIgCONrFnZrYfuKFBF9DtoEUCYWdJBdMWrGXpR3u44MxTuP+Ks0lum+B1WVKnKQIgH+h+2OtuQOGJXMDMXgVeTU1NvSWchYlI83ltXSG/emUDVTUh/uuHA5gwrLs+1eszTREAK4HezrmeQAEwHph4IhfQCEAkcu2rqObev27klZwCBnVP4sGrB9Hz5NZelyVH0dhtoPOA5UCqcy7fOXeTmdUAU4BMYBOQYWYbT+S6eiCMSGT64JO9XPynpSxeW8jU0b1ZcNsIdf4+1qgRgJlNqOf4EmBJY64tIpGjsqaWB97cwqx3PuG0Dq1YcNsIBvc4yeuy5Dh8eS8gTQGJRI4tu0q5c/4aNu3Yx4RhPbjnkr60TvBl1yJH8OV3SYvAIv4XChnPLMvj938/uL3zqWvTGF23vVMigy8DQET8bUdJOb/4y1re+3gvo/se3N55chtt74w0vgwATQGJ+Nfh2ztnjBvA+KHa3hmpfBkAmgIS8R9t7wweXwaAiPjL+5/s5ecZa9m5r4Kpo3szZVQvYmP0OJFI58sA0BSQiD9U1tTywBtbmLX0E1I6ttb2zoDxZQBoCkjEe7k7S5n60sHtnROHH9ze2Srel12GfEP6borIV4RCRmZeNS+/9S7tWsby9HVpXNBX2zuDSAEgIl/6/+2dVdreGQV8GQBaAxBpfq+uLeRXr6ynJmRc3y+ee69J0/bOgPPlMr5uBifSfErKq5k6P4c75uVwenIblvzkW6R3j1PnHwV8OQIQkeaxfOtefp6xhl2llfx0dB/+ddQZxMa0IM/rwqRZKABEolBlTS1/fGMLT9Zt73x58kgGdU/yuixpZgoAkSizeec+ps5fw+adpfxoeA9+pe2dUUvfdZEoEQoZs9/7lD/8PZd2ibHMvj6N756p7Z3RzJcBoF1AIuFVWHxwe+eyrXu58KxO3D9uAB21vTPq+TIA9ElgkfBZvLaQe+q2d/7+igFclaa7d8pBvgwAEWm8kvJq/v2vG/jrmkIG90jiwasGkaK7d8phFAAiAbRs6x5+kbGWXaWV/OzCPtyefobu3ilfowAQCZDDt3f27NiahZNHMlDbO6UeCgCRgDh8e+c15/bg7u9re6ccm/52iES4r27vjOOZ64cy6sxTvC5LIoAvA0DbQEUaRts7pTF8GQDaBipyfH9dU8CvF23Q9k75xnwZACJSv5Lyan69aAOL1xZyTo+DD2c/raO2d8qJUwCI+NiinAJmZuZSWFxOl6REfjC4C4tWF1BUWsnPL+zDZG3vlEZQAIj41KKcAqYvXE95dS0ABcXlPPbPrSS3iedlbe+UMNCPDiI+NTMz98vO/3CxMS3U+UtYKABEfKqwuPyox3eWVDRzJRJUCgARHyo+UEVC7NH/eXZJSmzmaiSoFAAiPvPBJ3u5+KGlVNWGiIv56rbOxLgYpo1J9agyCRoFgIhP1NSGePDNLUx48n1axsWweMr5zLxyIF2TEnFA16REZowbwNjBXb0uVQKi2XYBOedOB34FtDezK5vrfUUiQWFxOVPnr2FF3udccU43fvODfrRJiKV/1/bq8KXJNGgE4Jyb7Zwrcs5tOOL4Rc65XOfcx865u451DTP7xMxuakyxIkH09w07ufihpWwsLOFPVw/ij1cNpE2CdmhL02vo37I5wCPAs4cOOOdigEeBC4F8YKVzbjEQA8w44utvNLOiRlcrEiAV1bX89vUPef79bZzdrT0Pjx+sB7ZIs3Jm1rATnUsBXjOz/nWvRwD3mdmYutfTAczsyM7/yOssONYUkHNuEjAJIDk5eUhGRkaD6otEZWVltGnTxusymkSQ2waNb19BaYjH11aQX2Zc3DOOK3rHEdvCP/fx0fcvso0aNSrbzNKOd15jxpldge2Hvc4Hhtd3snOuI/A7YLBzbnp9QWFms4BZAKmpqZaent6IEv0tKyuLoLYvyG2Db94+M+PFFdv4j7c+pG3LOObeOIjv9EkOf4GNpO9fdGhMABztx5V6hxNmthe4rUEX1u2gJYBKDlRz18J1/G3DTr7V+2QeuGoQyW1162bxTmMCIB/oftjrbkBh48o5SLeDlqBZmfc5d87Loai0kru/fyY3n386LXw05SPRqTEBsBLo7ZzrCRQA44GJ4ShKIwAJitqQ8eg/P+ZPb22he4dWuomb+EpDt4HOA5YDqc65fOfcTWZWA0wBMoFNQIaZbQxHUWb2qplNCvIijQTfjpJyJj75Pg+8uYXLB3bhtTvOV+cvvtKgEYCZTajn+BJgSVgrEgmANz/cxbQFa6mqCfHHfxnIuHO66mld4ju+/LSJpoAkUlVU1zJjySbmLv+Mfl3a8T8TBnN6skay4k++DAAtAksk+riolCkv5rB5Zyk3n9+TaRelkhAb43VZIvXyZQCIRBIzI2PVdu5b/CGt4mN45vqhjDrzFK/LEjkuXwaApoAkUpSUV3P3K+t5fd0OzuvVkQevGsQp7Vp6XZZIg/gyADQFJJEg+7MvuHN+DjtLKvjlRWdy67e1t18iiy8DQMTPakPGq1urWPTGcjq3b0nGbSM4p8dJXpclcsJ8GQCaAhK/2rWvgp++tIZlW6u59OzO/Ne4AbRrGed1WSLfiC+fCKYPgokf/e/mXVz80FJythVzY/94/mfCYHX+EtF8OQIQ8ZPKmlru/9tmnnkvj76dD+7tz/9wlT7YJRFPASByDFt3l3HHizl8uGMf149M4a6Lz6RlXAz5H3pdmUjj+TIAtAYgXjMzFmTnc+/ijSTEtuCpa9MYfVYnr8sSCStfBoC2gYqXSiuq+dUrG1i8tpBzT+/An64ezKnttbdfgseXASDilTXbi7lj3moKiyv4xff6MDm9FzHa2y8BpQAQAUIhY9bST/jvzFw6tWtJxq3nMuS0Dl6XJdKkFAAS9YpKK/jZS2t59+M9fH/AqcwYdzbtE7W9U4JPASBRLSu3iJ9nrGV/VQ0zxg1g/NDu2t4pUcOXAaBdQNLUqmpCzMzczJNLP+XMU9syf8K59O7U1uuyRJqVLwNAu4CkKX26Zz8/mZfD+oISrh1xGnd/vy8t43Tffok+vgwAkaaycHU+v160gdiYFjzx4yGM6Xeq1yWJeEYBIFGhrLKGXy/awCs5BQzr2YE/XT2ILkmJXpcl4ikFgATeuvxifjIvh22fH+Cno/sw5bva2y8CCgAJsFDIePrdT/lD5maS2yQwf9IIhvXU3n6RQxQAEki7Syv5xV/W8vaW3Yzp14nfX3E2Sa3ivS5LxFd8GQDaBiqN8c6W3fwsYy2lFdX8dmx/fjS8h/b2ixyFLwNA20CloRblFDAzM5fC4nI6t2/JmZ3b8r+bd9P7lDa8cPNwUk/V3n6R+vgyAEQaYlFOAdMXrqe8uhaAwpIKCksqGHFGB2ZfN4zEeO3tFzkWXz4SUqQhZmbmftn5H27b3nJ1/iINoACQiFVYXH5Cx0XkqxQAEpHKq2ppGXf0v776gJdIwygAJOLsKCnnqieWU14dIvaID3QlxsUwbUyqR5WJRBYtAktEydn2BZOey+ZAZQ1PX5dGaUXNl7uAuiQlMm1MKmMHd/W6TJGIoACQiPFKTj6/fHk9p7ZryQs3D6dP3e2b1eGLfDMKAPG9UMiY+UYuj2dtZXjPDjx+zRA6tNanekUaq9kCwDk3FrgEOAV41MzeaK73lshVVlnD1Pk5vLWpiInDe3DfZf2Ij9XSlUg4NOhfknNutnOuyDm34YjjFznncp1zHzvn7jrWNcxskZndAlwPXP2NK5aosf3zA1zx2DL+mbub//hBP343tr86f5EwaugIYA7wCPDsoQPOuRjgUeBCIB9Y6ZxbDMQAM474+hvNrKju1/fUfZ1Ivd7/ZC+Tn8+mNmTMvWEY5/c+2euSRALHmVnDTnQuBXjNzPrXvR4B3GdmY+peTwcwsyM7/0Nf74D7gTfN7K1jvM8kYBJAcnLykIyMjAY3JtKUlZXRpk0br8toEo1pW9b2ap77sIrkVo6p57Tk1Nb++6k/yN87UPsi3ahRo7LNLO145zVmDaArsP2w1/nA8GOcfwcwGmjvnOtlZn8+2klmNguYBZCammrp6emNKNHfsrKyCGr7vknbampD/Pb1TczZmMd3+iTz8ITBtE+Ma5oCGynI3ztQ+6JFYwLgaPfXrXc4YWYPAw836MK6HXTUKTlQzZR5q1n60R5uPr8n07/fV0/tEmlijRlb5wPdD3vdDShsXDkHmdmrZjYpyEM0+X9bd5cx9rH3eP+TvfzhirO559Kz1PmLNIPGjABWAr2dcz2BAmA8MDEcRWkEED3e3rKbKS+uJj6mBS/eci5DU/TIRpHm0tBtoPOA5UCqcy7fOXeTmdUAU4BMYBOQYWYbw1GURgDBZ2bMfvdTbnhmBV2TEvnrlPPU+Ys0swaNAMxsQj3HlwBLwlqRBF5VTYhfL9rAS6u2872zOvHg1YNonaAPpYs0N1/+q9MUUHDtLatk8vOrWZH3OXd8txc/Hd2HFprvF/GE/zZYoymgoNq0Yx+XP/Iea/OLeWj8IH7+vVR1/iIe8uUIQILnjY07mfrSGtq2jCXj1hEM7J7kdUkiUc+XAaApoOAwMx7L2srMzFwGdmvPrGvT6NSupddliQiaApImVFFdy9SX1jAzM5fLB3bhpVtHqPMX8RFfjgAk8n1REeLqJ5azNr+EaWNSuT39DA7eDkpE/MKXAaApoMi2dnsxv1leQZVVMevHQ/hev1O9LklEjkJTQBJWf11TwFVPLCfGwcuTR6rzF/ExX44AJPKEQsYDb27hkX9+zLCUDlxzegV9O7fzuiwROQZfjgAksuyvrOG257N55J8fc3Vad56/eTjt4jXfL+J3vhwBaA0gcuR/cYCb565iy65S/v3Ss7jhvBQt9opECF+OALQGEBlW5n3ODx55j4Licp65YRg3nt9Tnb9IBPHlCED876WV27hn0Qa6ndSKp65L44xkhbVIpFEAyAmpqQ0x42+befrdTzm/18k8OvEc2rfy52MbReTYFADSYCXl1dwxL4d3tuzm+pEp3HNJX2JjfDmLKCINoACQBvl0z35umruSbXsP8F8/HMDE4T28LklEGsmXAaBdQP7y7kd7uP2FbGJaOJ6/eTjnnt7R65JEJAx8OX7XLiB/MDPmLsvjumdWcGr7liyecr46f5EA8eUIQLxXVRPi3sUbmbdiG6P7nsKDVw+ibUst9ooEiQJAvubz/VVMfj6bDz79nMnpZ/CL76USoyd3iQSOAkC+YsuuUm6au5Jd+yp58OqB/HBwN69LEpEmogCQL7314S7unJ9Dq4RYXpp0LoN7nOR1SSLShBQAgpnxxDuf8Pu/b6Zfl3Y8eW0andsnel2WiDQxXwaAtoE2n4rqWu5euJ6FOQVccnZn/vvKgSTGx3hdlog0A20DjWJF+yoYP+t9FuYU8LML+/DIhMHq/EWiiC9HANL0NhSUcMuzqyg+UM3jPzqHiwd09rokEWlmCoAo9Pq6Hfz8L2vo0CqeBZNH0K9Le69LEhEPKACiwKKcAmZm5lJYXE6blrGUVtRwTo8knvhxGsltE7wuT0Q8ogAIuEU5BUxfuJ7y6loASitqiHGOicN6qPMXiXK+XASW8JmZufnLzv+QWjMefOsjjyoSEb/QCCCgyiprWLg6n4LiiqP+fmFxeTNXJCJ+owAImE/37Gfusjxezs6ntLKGuBhHda197bwuSfqgl0i0UwAEQChkvP3RbuYuyyMrdzdxMY7vD+jM9SNT+Gzvga+sAQAkxsUwbUyqhxWLiB80WwA45/oCdwInA/8ws8eb672Dal9FNQtW5fPc+5/x6Z79JLdNYOro3kwc3oNT2rYE+PJ+Pod2AXVJSmTamFTGDu7qZeki4gMNCgDn3GzgUqDIzPofdvwi4CEgBnjKzO6v7xpmtgm4zTnXAniyUVVHuY+Lynh2+cFpnv1VtQzukcRD4wdxcf/OxMd+fV1/7OCu6vBF5GsaOgKYAzwCPHvogHMuBngUuBDIB1Y65xZzMAxmHPH1N5pZkXPucuCuumvJCagNGVm5RcxZlsfSj/YQH9OCSwcenOY5u1uS1+WJSARyZl9fIDzqic6lAK8dGgE450YA95nZmLrX0wHM7MjO/2jXet3MLqnn9yYBkwCSk5OHZGRkNKi+SFRWVsbx7ne0v9pYml/DP7ZVs7vcSEpwfLdHLOnd4miX4N+HtDSkbZFM7YtsQW/fqFGjss0s7XjnNWYNoCuw/bDX+cDw+k52zqUD44AEYEl955nZLGAWQGpqqqWnpzeiRH/LysqivvZt2VXKnGV5vLK6gPLqWoamnMS9I1MY0+9U4mL8//GNY7UtCNS+yBb09jVUYwLgaD9+1jucMLMsIKtBF47S20HXhoy3Nu1i7rI8lm3dS3xsC8YO6sK1I1Lo31X36xGR8GpMAOQD3Q973Q0obFw5B5nZq8Crqampt4Tjen5XfKCK+Su389zyzygoLqdL+5b820WpjB/agw6t470uT0QCqjEBsBLo7ZzrCRQA44GJYakqSmwvDXHXy+t4JaeAypoQ557egV9f2pfRfTsRGwHTPCIS2Rq6DXQekA6c7JzLB+41s6edc1OATA7u/JltZhvDUVSQp4BqakO88eEu5izLY8Wn5bSMK2DcOV25dkQKfTu387o8EYkiDQoAM5tQz/ElHGNB95sK4hTQ3rJK5q/czvPvf8aOkgq6nZTI1anxTL/6OyS10jSPiDQ/X94KIkgjgA0FJcxZlsfitYVU1YQ4r1dHfnN5Py7o24ml77ytzl9EPOPLAIj0EUB1bYi/b9jJnGV5ZH/2Ba3iY7gqrRvXjUihd6e2XpcnIgL4NAAi1e7SSuat2MYLH3zGrn2VnNaxFfdc0pd/SetO+8Q4r8sTEfkKXwZApE0BrdlezNxleby+bgdVtSG+3SeZGeNOI73PKbRo4d9P64pIdPNlAETCFFBVTYgl63fwzLI81m4vpnV8DBOGdefakSmckRzcj5iLSHD4MgD8rGhfBc9/sI0XP9jGnrJKTj+5NfdddhZXDOlG25aa5hGRyOHLAPDbFJCZsXrbwWmeJet3UBMyRqUmc93IFL7dO1nTPCISkXwZAH6ZAqqoruW1dTuYuyyP9QUltE2I5doRKVw74jRSTm7tZWkiIo3mywBobotyCr7yxKxbvt2TPaVVzFuxjb37q+h1Shv+c2x/xg3uSusE/ZGJSDBEfW+2KKfgK8/MLSgu577FHwIwum8nrh+Zwnm9OuKcpnlEJFiiPgBmZuZ+5YHph3Rql8BT1x33eQoiIhHLl7ecdM5d5pybVVZW1uTvVVhcftTjRfsqm/y9RUS85MsAMLNXzWxSczyyrUtS4gkdFxEJCl8GQHOaNiaVxLiYrxxLjIth2phUjyoSEWkeUb8GMHZwV4Cv7AKaNib1y+MiIkEV9QEAB0NAHb6IRBtfTgE15yKwiEi08mUANOcisIhItPJlAIiISNNTAIiIRCkFgIhIlHJm5nUN9XLOlQK5XtfRhE4G9nhdRBMJcttA7Yt0QW9fqpkd9wHkft8Gmmtmgb0hj3NuVVDbF+S2gdoX6aKhfQ05T1NAIiJRSgEgIhKl/B4As7wuoIkFuX1BbhuofZFO7cPni8AiItJ0/D4CEBGRJqIAEBGJUr4MAOfcbOdckXNug9e1hJtzrrtz7p/OuU3OuY3OuTu9rimcnHMtnXMrnHNr69r3G69rCjfnXIxzLsc595rXtTQF51yec269c25NQ7cTRgrnXJJzboFzbnPdv8ERXtcULs651Lrv2aH/9jnnph7za/y4BuCc+zZQBjxrZv29riecnHOdgc5mtto51xbIBsaa2YcelxYWzjkHtDazMudcHPAucKeZve9xaWHjnPsZkAa0M7NLva4n3JxzeUCamQXug1LOubnAUjN7yjkXD7Qys2KPywo751wMUAAMN7PP6jvPlyMAM3sH+NzrOpqCme0ws9V1vy4FNgGBeRiBHXToPt5xdf/576eMb8g51w24BHjK61rkxDjn2gHfBp4GMLOqIHb+dS4Ath6r8wefBkC0cM6lAIOBDzwuJazqpkjWAEXAm2YWpPb9Cfg3IORxHU3JgDecc9nOuUleFxNGpwO7gWfqpvCecs619rqoJjIemHe8kxQAHnHOtQFeBqaa2T6v6wknM6s1s0FAN2CYcy4Q03jOuUuBIjPL9rqWJnaemZ0DXAz8a92UbBDEAucAj5vZYGA/cJe3JYVf3dTW5cBfjneuAsADdXPjLwMvmNlCr+tpKnXD6yzgIm8rCZvzgMvr5sjnA991zj3vbUnhZ2aFdf8vAl4BhnlbUdjkA/mHjUgXcDAQguZiYLWZ7TreiQqAZla3SPo0sMnMHvC6nnBzziU755Lqfp0IjAY2e1pUmJjZdDPrZmYpHBxi/6+ZXeNxWWHlnGtdtzmBuumR7wGB2I1nZjuB7c651LpDFwCB2HxxhAk0YPoHfHo3UOfcPCAdONk5lw/ca2ZPe1tV2JwH/BhYXzdPDnC3mS3xrqSw6gzMrduF0ALIMLNAbpcMqE7AKwd/TiEWeNHM/u5tSWF1B/BC3TTJJ8ANHtcTVs65VsCFwK0NOt+P20BFRKTpaQpIRCRKKQBERKKUAkBEJEopAEREopQCQEQkSikARESilAJARCRK/R/NjVDn8BgVdQAAAABJRU5ErkJggg==\n", "text/plain": [ "
\n", "
" ] }, { "cell_type": "markdown", "id": "03709f67-97c5-4cc0-8c80-76bda8027a92", "metadata": {}, "source": [ "**Question 2.5** Not a very high probability, perhaps, but then a year isn’t very long! Let’s say you are living in the Los Angeles area and take out a 30-year mortgage on a house. What is the probability of a magnitude 7.0-7.9 earthquake occurring during that 30-year period? How does this value compare to the probability of such an earthquake occurring in the San Francisco area in the next 30 years? *Hint:* You can modify the code for calculating the joint probability for San Francisco as provided below:" ] }, { "cell_type": "code", "execution_count": 86, "id": "4faf8bd9-3ac7-44b1-9fad-99719c0466a0", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Probability of a magnitude 7.0-7.9 earthquake occurring in the San Francisco are during that 30-year period is 45.45 percent.\n" ] } ], "source": [ "MRI = 50\n", "num_years = 30\n", "prob = (1 - (1 - 1/MRI) ** num_years)*100\n", "print(\"Probability of a magnitude 7.0-7.9 earthquake occurring in the San Francisco are during that 30-year period is %.2f percent.\"%(prob)) " ] }, { "cell_type": "markdown", "id": "93304e6a-5c7f-4a40-86c3-c2df7b6110ac", "metadata": {}, "source": [ "**Answer:**" ] }, { "cell_type": "markdown", "id": "bf5d1df8-ec63-46c9-8a24-51089e7bc950", "metadata": {}, "source": [ "
\n", "
" ] }, { "cell_type": "markdown", "id": "cb4fde1c-582e-4309-b6f4-78740c68bffa", "metadata": {}, "source": [ "----\n", "### Mini-tutorial\n", "\n", " - **Mean and Standard Deviation with numpy**\n", " - **Statistical distributions and confidence intervals**\n", " \n", "Consider a set of observations $x_i$ from $i = 1$ to $N$. The mean observation is denoted by $\\bar{x}$ and is obtained by summing all observations $x_i$ and dividing by the number of observation $N$:\n", "\n", "$$ \\bar{x} = \\frac{1}{N}\\sum_{i=1}^N x_i = \\frac{x_1 + x_2 + x_3 + \\cdots + x_N}{N} $$Standard deviation ($\\sigma$) is a measure that represents about the spread in observed values in the sample of observations. The formula is not intuitive to most people:\n", "\n", "$$ \\sigma = \\sqrt{\\frac{\\sum_{i=1}^N (x_i - \\bar{x})^2}{N - 1}} $$\n", "\n", "The equation can be read as: *For each of the $N$ observations, calculate the deviation of the\n", "observation from the calculated mean. Square all of these deviations and sum them. Divide this sum by the total number of observations $N$ minus one. Finally, take the\n", "square root of this number to arrive at the standard deviation $\\sigma$.*" ] }, { "cell_type": "code", "execution_count": null, "id": "c8cc2a6e-1c86-48f8-aad2-4d37f2ecc3b1", "metadata": {}, "outputs": [], "source": [ "from scipy import stats, optimize, special, integrate\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": null, "id": "a6fffc51-8897-4111-abed-bbac93e5e6e7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x = [ 2 3 5 10 15]\n", "mean = 7.00\n", "std = 4.86\n" ] } ], "source": [ "x = np.array([2, 3, 5, 10, 15])\n", "print(\"x =\", x)\n", "print(\"mean = {:.2f}\".format(np.mean(x)))\n", "print(\"std = {:.2f}\".format(np.std(x)))" ] }, { "cell_type": "markdown", "id": "fb57f33d-7425-436c-ae67-fe50f235c167", "metadata": {}, "source": [ "Once you know the mean $\\bar{x}$ and standard deviation $\\sigma$ of a *normal distribution*, you can\n", "define the probability distribution that describes the probability of obtaining a measurement\n", "in a certain range of $x$. Note that all probability distributions have unit area such that $\\int_{-\\infty}^{\\infty} P(x)$ = 1. The probability distribution is a curve that can be described\n", "mathematically by the formula\n", "\n", "$$ P(x) = \\frac{1}{\\sigma \\sqrt{2\\pi}} \\exp \\left\\{-\\frac{(x - \\bar{x})^2}{2\\sigma^2} \\right\\} $$\n", "\n", "\n", "In this formula, $P(x)$ does not represent the probability of observing exactly $x$. Instead, $P(x)$\n", "is called a probability density function, and the probability of observing a value between two\n", "arbitrary values $a$ and $b$ is the area beneath the curve in the range $[a, b]$ (meaning 'a to b' with both values included).\n", "For example, the probability that the next observation x lies between the mean minus one\n", "standard deviation and the mean plus one standard deviation (between $\\bar{x}-\\sigma$ and $\\bar{x}+\\sigma)$ is \n", "simply the area beneath $P(x)$ over the interval $[\\bar{x}-\\sigma, \\bar{x}+\\sigma]$. This area turns out to be 0.68.\n", "Since 0.68 = 68%, you can also express this result as: *we are 68% confident that the measurement will lie between $\\bar{x}-\\sigma$ and $\\bar{x}+\\sigma$*.\n", "\n", "It is important to remember that quite often, we do not actually observe the normal distribution. We\n", "collect a *sample* (usually not that large) and calculate the mean and the standard deviation.\n", "We then make the assumption that the *underlying population* has a normal distribution. This\n", "is not always true. When only a small number of measurements from a *sample* are used to estimate $\\bar{x}$ and $\\sigma$ of the entire *population*, these numbers are less certain. In such cases, the probability distribution is \n", "more accurately described by *Student's t-distribution*, which takes into account the\n", "uncertainties in $\\bar{x}$ and $\\sigma$ as well due to the limited sampling. Basically, this uncertainty will act to *spread out* the\n", "normal distribution when only few observations are used, since we are less certain about the distribution." ] }, { "cell_type": "markdown", "id": "456e5ef7-674e-40a8-b74a-baf043ade699", "metadata": { "tags": [] }, "source": [ "
![](\"Images/PDF_comparison.png\")
"
]
},
{
"cell_type": "markdown",
"id": "f3b1cbd4-b55a-441c-87d1-a78053c161d7",
"metadata": {},
"source": [
"The probability distribution of a *Student's t-distribution* is a curve that has unit area and is described mathematically by the formula:\n",
"\n",
"$$ P(x|\\nu, \\bar{x}, \\sigma) = \\frac{\\Gamma\\left(\\frac{\\nu+1}{2}\\right)}{\\sigma\\sqrt{\\nu\\pi}\\ \\Gamma\\left(\\frac{\\nu}{2}\\right)} \\left(1 + \\frac{(x-\\bar{x})^2}{\\sigma^2 \\nu} \\right)^{-\\frac{\\nu+1}{2}} $$\n",
"\n",
"where $\\nu$ is the degree of freedom (number of samples minus one), $\\bar{x}$ is the sample mean, $\\sigma$ is the sample standard deviation, and $\\Gamma(z) = (z - 1)!$ (see [the gamma function](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.gamma.html))."
]
},
{
"cell_type": "markdown",
"id": "fab4cedc-7480-4999-885c-aa4f94d57d51",
"metadata": {},
"source": [
"The example below illustrates how samples are generated from a Normal distribution with known population mean and population standard deviation. When you answer the questions below, you will be supplied with the data (treated as *samples*) without any knowledge of the underlying *population* mean and population standard deviation."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b219b3e2-1bba-4c66-b901-9baa55f192f3",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([49.30996753, 20.97160251, 40.41776131, 18.44689402, 44.1202185 ])"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"popmean,popstd = 30, 12\n",
"N = 5\n",
"#s = np.random.normal(popmean, popstd)\n",
"s = np.array([49.30996753, 20.97160251, 40.41776131, 18.44689402, 44.1202185 ]) # one random output from the command above\n",
"s"
]
},
{
"cell_type": "markdown",
"id": "ba67809b-3fb4-45a3-a8aa-6b45cf992723",
"metadata": {},
"source": [
"Let's pretend we do not know the population mean and the population standard deviation, and then we calculate the sample mean and the sample standard deviation."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d8a0069a-b762-43e0-b55f-77bc83d22ef8",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Sample size N = 5\n",
"Sample mean = 34.65\n",
"Sample std = 12.55\n"
]
}
],
"source": [
"samplemean = np.mean(s)\n",
"samplestd = np.std(s)\n",
"\n",
"print(\"Sample size N = %d\" % N)\n",
"print(\"Sample mean = %.2f\" % samplemean)\n",
"print(\"Sample std = %.2f\" % samplestd)"
]
},
{
"cell_type": "markdown",
"id": "b925a8ff-e0fc-4604-a6e9-e1ef3a131407",
"metadata": {},
"source": [
"Now, we will construct the probability density of the Student's t-distribution from the *sample mean, the *sample* standard deviation, and the *sample* size using the function below."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5eed5ae7-fef2-40dc-b929-0108c806c50e",
"metadata": {},
"outputs": [],
"source": [
"def tdist_probability_density(x, dof, mean, std):\n",
" '''\n",
" Compute the probability density of the Student's t-distribution\n",
" \n",
" Input Parameters:\n",
" ----------------\n",
" \n",
" x: (array of) location(s) where the probability density is computed\n",
" dof: degree of freedom; equals to the sample size minus one\n",
" mean: sample mean\n",
" std: sample standard deviation\n",
" \n",
" Return:\n",
" ------\n",
" \n",
" p: probability density at all given x\n",
" \n",
" '''\n",
" p = special.gamma((dof + 1) / 2) / special.gamma(dof / 2) / np.sqrt(np.pi * dof * std ** 2) * \\\n",
" (1 + (x - mean) ** 2 / std ** 2 / dof) ** (-(dof + 1) / 2)\n",
" return p\n",
"\n",
"def ndist_probability_density(x, mean, std):\n",
" '''\n",
" Compute the probability density of the Normal distribution\n",
" \n",
" Input Parameters:\n",
" ----------------\n",
" \n",
" x: (array of) location(s) where the probability density is computed\n",
" mean: sample mean\n",
" std: sample standard deviation\n",
" \n",
" Return:\n",
" ------\n",
" \n",
" p: probability density at all given x\n",
" \n",
" '''\n",
" p = np.exp(-1/2 * (x - mean) ** 2 / std ** 2) / np.sqrt(2 * np.pi * std ** 2)\n",
" return p"
]
},
{
"cell_type": "markdown",
"id": "bfd5efc2-055a-428e-94f4-5f62e131efc3",
"metadata": {},
"source": [
"The *cumulative density function* (CDF) is obtained by integrating the probability density function from $-\\infty$ to $x$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a33f072d-c518-4b07-abdc-5ea56a5a376f",
"metadata": {},
"outputs": [],
"source": [
"def tdist_cumulative_density(x, dof, mean, std):\n",
" '''\n",
" Compute the cumulative density of the Student's t-distribution\n",
" \n",
" Input Parameters:\n",
" ----------------\n",
" \n",
" x: (array of) location(s) where the probability density is computed\n",
" dof: degree of freedom; equals to the sample size minus one\n",
" mean: sample mean\n",
" std: sample standard deviation\n",
" \n",
" Return:\n",
" ------\n",
" \n",
" c: cumulative density at all given x\n",
" \n",
" '''\n",
" c = integrate.quad(tdist_probability_density, -np.inf, x, args=(dof, mean, std))[0]\n",
" return c\n",
"\n",
"def ndist_cumulative_density(x, mean, std):\n",
" '''\n",
" Compute the cumulative density of the Normal distribution\n",
" \n",
" Input Parameters:\n",
" ----------------\n",
" \n",
" x: (array of) location(s) where the probability density is computed\n",
" mean: sample mean\n",
" std: sample standard deviation\n",
" \n",
" Return:\n",
" ------\n",
" \n",
" c: cumulative density at all given x\n",
" \n",
" '''\n",
" c = integrate.quad(ndist_probability_density, -np.inf, x, args=(mean, std))[0]\n",
" return c"
]
},
{
"cell_type": "markdown",
"id": "5cca4cde-697b-48fe-9e3b-3434351d8b2b",
"metadata": {},
"source": [
"Here is how to calculate the value such that there is a 95% chance that a sample drawn from the Student's t-distribution is smaller"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0caeba1f-e8aa-4298-bbe6-c6802c3428e7",
"metadata": {},
"outputs": [],
"source": [
"def find_confidence_interval_onesided(dof, mean, std, prob, distribution='t'):\n",
" '''\n",
" Computes the value to give any probability that a sample draw from \n",
" the Student's t-distribution is smaller than this value. In other words,\n",
" this function finds the value such that CDF(value | dof, mean, std) = prob.\n",
" \n",
" Input Parameters:\n",
" ----------------\n",
" \n",
" dof: degree of freedom; equals to the sample size minus one\n",
" mean: sample mean\n",
" std: sample standard deviation\n",
" prob: target probability\n",
" distribution: which type of distribution ('t' or 'normal')\n",
" \n",
" Return:\n",
" ------\n",
" \n",
" value: a value such that CDF(value) = prob\n",
" \n",
" '''\n",
" # function to optimize\n",
" # which is the difference between the cumulative distribution and target probability\n",
" if distribution == 't':\n",
" def func_to_optimize(x, dof, mean, std, prob):\n",
" return tdist_cumulative_density(x, dof, mean, std) - prob\n",
" else:\n",
" def func_to_optimize(x, dof, mean, std, prob):\n",
" return ndist_cumulative_density(x, mean, std) - prob\n",
" \n",
" # x0 and x1 are initial and second guesses\n",
" value = optimize.root_scalar(func_to_optimize, args=(dof, mean, std, prob), x0=mean, x1=mean+std).root\n",
" return value"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ca539a44-86f8-4851-9bbe-e0534033593e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"61.407853899787334"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"find_confidence_interval_onesided(N-1, samplemean, samplestd, 0.95) # 95% == 0.95"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6bb06350-c36e-47bd-af50-907887827ca9",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4.604094871350386"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P = 0.005\n",
"find_confidence_interval_onesided(N-1, 0, 1, 1-P) "
]
},
{
"cell_type": "markdown",
"id": "1792604a-e320-4a8d-a242-c7d0d8d4d88c",
"metadata": {},
"source": [
"Values of $t$ score ($t = \\frac{x - \\bar{x}}{\\sigma}$) corresponding to given values of $P$ and $\\nu$ degree of freedom are provided in the table below. The value $P$ is defined as the area under the probability density curve to the right of $t$. The area under the probability density curve to the left is the cumulative density function. The value corresponding to the figure below is marked in **bold** inside the table.\n",
"\n",
"![](\"Images/CDF_and_P.png\")
\n",
"\n",
"| Degree of freedom $\\nu$ |P=0.005|P=0.01 |P=0.025|P=0.05 |P=0.10 |P=0.25 |\n",
"|-------------------------|-------|-------|-------|-------|-------|-------|\n",
"| 1 |63.657 |31.281 |12.706 | 6.314 | 3.078 | 1.000 |\n",
"| 2 | 9.925 | 6.945 | 4.303 | 2.920 | 1.886 | 0.816 |\n",
"| 3 | 5.841 | 4.541 | 3.182 | 2.353 | 1.638 | 0.765 |\n",
"| 4 | 4.604 | 3.747 | 2.776 | 2.132 | **1.533** | 0.741 |\n",
"| 5 | 4.032 | 3.365 | 2.571 | 2.015 | 1.476 | 0.727 |\n",
"| 6 | 3.707 | 3.143 | 2.447 | 1.943 | 1.440 | 0.718 |\n",
"| 7 | 3.499 | 2.998 | 2.365 | 1.895 | 1.415 | 0.711 |\n",
"| 8 | 3.355 | 2.896 | 2.306 | 1.860 | 1.397 | 0.706 |\n",
"| 9 | 3.250 | 2.821 | 2.262 | 1.833 | 1.383 | 0.703 |\n",
"| 10 | 3.169 | 2.764 | 2.228 | 1.812 | 1.372 | 0.700 |\n",
"| 15 | 2.947 | 2.602 | 2.131 | 1.753 | 1.341 | 0.691 |\n",
"| 20 | 2.845 | 2.528 | 2.086 | 1.725 | 1.325 | 0.687 |\n",
"| 25 | 2.787 | 2.485 | 2.060 | 1.708 | 1.316 | 0.684 |\n",
"| 30 | 2.750 | 2.457 | 2.042 | 1.697 | 1.310 | 0.683 |\n",
"| $\\infty$ | 2.576 | 2.326 | 1.960 | 1.645 | 1.282 | 0.674 |"
]
},
{
"cell_type": "markdown",
"id": "ae6388de-e778-4290-b238-cb78222b4969",
"metadata": {
"tags": []
},
"source": [
"----\n",
"\n",
"## Part III: The Parkfield Prediction\n",
"\n",
"It is important to realize that the general observation that earthquakes recur on individual faults can be the result of many physical realities. The simple ‘sawtooth’ earthquake cycle model may be attractive because of its simplicity, but its validity has to be confirmed by observational data. This model has come under a great deal of criticism lately, because its predictions for the locations of earthquakes over the past 20 years along the Pacific rim have not been very successful. "
]
},
{
"cell_type": "markdown",
"id": "26b313cd-8072-4059-aa11-03bea6920492",
"metadata": {},
"source": [
"![](\"Images/Parkfield_map.png\")
"
]
},
{
"cell_type": "markdown",
"id": "a2dfd3a9-849a-4e80-b9ec-53b06d3a7a6c",
"metadata": {
"tags": []
},
"source": [
"The purpose of this problem is to illustrate how probabilistic forecasts are made by calculating the likely time for next Parkfield earthquake using statistics.\n",
"\n",
"*Background*: The National Earthquake Prediction Evaluation Council (NEPEC) in 1985 endorsed a prediction that with 95% confidence (= probability), the next Parkfield earthquake would occur\n",
"before January 1993."
]
},
{
"cell_type": "markdown",
"id": "42104e46-25f6-4faa-aa1f-5a1f56562c1a",
"metadata": {},
"source": [
"![](\"Images/Parkfield_events.png\")
"
]
},
{
"cell_type": "markdown",
"id": "a0585467-276c-4ae2-a7cd-35332dd137f4",
"metadata": {
"tags": []
},
"source": [
"*The Data*: At the time, most seismologists agreed that earthquakes of very similar size and character (that is, which ruptured the same part of the San Andreas Fault) have occurred near Parkfield at the following times: 1857, 1881, 1901, 1922, 1934, and 1966. There is no record of earthquake activity before 1857, but we do not expect there to be one because the area was not settled at that time.\n",
"\n",
"In order to calculate expectations and probabilities, you have to have a model in mind which tells you what the data above represent. It is not clear which model is most appropriate, and you may think of a better one than the following:\n",
"\n",
"*Approach 1 — strictly numbers.* The sequence clearly suggests a repetition with approximately the same interval between events. If we think of the intervals between earthquakes as belonging to a normal distribution with some mean, the data can be used to estimate both an average repeat time and a standard deviation around that mean."
]
},
{
"cell_type": "markdown",
"id": "6803ef1f-ead5-41ee-b99d-6d9d2a6cadae",
"metadata": {},
"source": [
"**Question 3.1** Calculate the mean value and the standard deviation for the earthquake repeat times. *Hint:* Create a numpy array for `years = np.array([1857, 1881, 1901, 1922, 1934, 1966])` and use standard numpy functions."
]
},
{
"cell_type": "markdown",
"id": "74cdd789-e760-462a-87dc-5d6704770243",
"metadata": {},
"source": [
"**Answer:**"
]
},
{
"cell_type": "markdown",
"id": "e0780d76-d27e-4522-879a-5834ee20e8f8",
"metadata": {},
"source": [
"\n", "
" ] }, { "cell_type": "markdown", "id": "3cd0f9ec-5c28-4d75-941a-1901006ff3b5", "metadata": {}, "source": [ "**Question 3.2** Based on the values derived above, what is the predicted (i.e. most likely) time for an earthquake to occur following the 1966 earthquake ?" ] }, { "cell_type": "markdown", "id": "f81b94bf-e673-4231-a5b4-e731fbb7e079", "metadata": {}, "source": [ "**Answer:**" ] }, { "cell_type": "markdown", "id": "938b4007-28cf-48f8-a540-79fd6dbe8ef8", "metadata": {}, "source": [ "
\n", "
" ] }, { "cell_type": "markdown", "id": "e1110b28-8a0d-4c44-b1a0-5bf13198e489", "metadata": {}, "source": [ "**Question 3.3** What is the time before which there is a 95% probability that the earthquake will occur? *Hint:* Since there are very few observations, you need to know something about Student’s t-distributions to answer this question. See **mini-tutorial** above." ] }, { "cell_type": "markdown", "id": "a27f9f51-cbc7-4745-912d-567acbcd8ee1", "metadata": {}, "source": [ "**Answer:**" ] }, { "cell_type": "markdown", "id": "7db04290-2035-4d69-b306-08efff57fc23", "metadata": {}, "source": [ "
\n", "
" ] }, { "cell_type": "markdown", "id": "d87c0b61-c9fd-460f-a06b-c86df561d092", "metadata": {}, "source": [ "*Approach 2 — using some ‘insights’.* Assume that the 1934 earthquake really should have occurred in 1944, but that it for some reason was ‘triggered’ to occur prematurely by some other geophysical process." ] }, { "cell_type": "markdown", "id": "bdd50d83-613e-437f-9432-afaa1b13a44a", "metadata": {}, "source": [ "**Question 3.4** Repeat the analysis for mean and standard deviation but assume that the 1934 event really occurred in 1944. " ] }, { "cell_type": "markdown", "id": "abf1da3d-67b5-4274-b869-42ddfb507599", "metadata": {}, "source": [ "**Answer:**" ] }, { "cell_type": "markdown", "id": "da2a1887-3b60-4b9f-82d7-7a8e28df4960", "metadata": {}, "source": [ "
\n", "
" ] }, { "cell_type": "markdown", "id": "5f1c9cf2-0f7a-49da-a604-8258ff53561e", "metadata": {}, "source": [ "**Question 3.5** Which approach appears to be closest to the one used by the NEPEC? If we were uncertain about the correctness of different approaches to interpreting the data (for example Approaches 1 and 2 above) how might one incorporate these uncertainties in the probability estimates?" ] }, { "cell_type": "markdown", "id": "febfa9e3-d60f-486b-95a0-97612ee69479", "metadata": {}, "source": [ "**Answer:**" ] }, { "cell_type": "markdown", "id": "2e1e7b2c-1e2a-4c7f-9e96-413e58de2f77", "metadata": {}, "source": [ "
\n", "
" ] }, { "cell_type": "markdown", "id": "9d8c8bac-ea78-4ba3-9d12-72f92da4609a", "metadata": {}, "source": [ "**Question 3.6** The long-awaited Parkfield earthquake finally occurred on September 28, 2004! Include the 2004 data point to \n", "\n", "(a) Re-calculate the mean value and standard deviation for the earthquake repeat time with this year.\n", "\n", "(b) Predict the year for the next Parkfield earthquake with 95 % probability." ] }, { "cell_type": "markdown", "id": "9fde297e-64d4-49a6-9924-63285c915104", "metadata": {}, "source": [ "**Answer:**" ] }, { "cell_type": "markdown", "id": "db485bd4-1d29-45c6-aa4f-d09d4559fb6a", "metadata": {}, "source": [ "
\n", "
" ] } ], "metadata": { "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.9.12" } }, "nbformat": 4, "nbformat_minor": 5 }