From 8a1dc5ea34d7b9ab52a257671304231e26f309c7 Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Mon, 16 Apr 2018 07:01:31 -0400 Subject: [PATCH] 2.11.13 in progress --- BDA 2.11.13.ipynb | 438 ++++++++++++++++++++++++++++++++++++++++++++++ BDA 2.11.13.md | 32 ++++ BDA 2.11.13.txt | 9 - 3 files changed, 470 insertions(+), 9 deletions(-) create mode 100644 BDA 2.11.13.ipynb create mode 100644 BDA 2.11.13.md delete mode 100644 BDA 2.11.13.txt diff --git a/BDA 2.11.13.ipynb b/BDA 2.11.13.ipynb new file mode 100644 index 0000000..46867cc --- /dev/null +++ b/BDA 2.11.13.ipynb @@ -0,0 +1,438 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Problem 2.10.13\n", + "\n", + "Discrete data: Table 2.2 gives the number of fatal accidents and deaths on scheduled airline flights per year over a ten-year period. \n", + "We use these data as a numerical example for fitting discrete data models. \n", + "\n", + "1. Assume that the numbers of fatal accidents in each year are independent with a Poisson(theta) distribution. Set a prior distribution for theta and determine the posterior distribution based on the data from 1976 through 1985. Under this model, give a 95% predictive interval for the number of fatal accidents in 1986. You can use the normal approximation to the gamma and Poisson or compute using simulation.\n", + "2. Assume that the numbers of fatal accidents in each year follow independent Poisson distributions with a constant rate and an exposure in each year proportional to the number of passenger miles flown. Set a prior distribution for theta and determine the posterior distribution based on the data for 1976–1985. (Estimate the number of passenger miles flown in each year by dividing the appropriate columns of Table 2.2 and ignoring round-off errors.) Give a 95% predictive interval for the number of fatal iaccidents in 1986 under the assumption that 8 × 10 11 passenger miles are flown that year.\n", + "3. Repeat (1) above, replacing ‘fatal accidents’ with ‘passenger deaths.’\n", + "4. Repeat (2) above, replacing ‘fatal accidents’ with ‘passenger deaths.’\n", + "5. In which of the cases above does the Poisson model seem more or less reasonable? Why? Discuss based on general principles,without specific reference to the numbers in Table 2.2. Incidentally, in 1986, there were 22 fatal accidents, 546 passenger deaths, and a death rate of 0.06 per 100 million miles flown. We return to this example in Exercises 3.12, 6.2, 6.3, and 8.14.\n", + "\n", + "|Year |Fatal accidents |Passenger deaths |Death rate\n", + "|---|---|---|---| \n", + "|1976 | 24 | 734 | 0.19 \n", + "|1977 |25 |516 |0.12 \n", + "|1978 |31 |754 |0.15 \n", + "|1979 |31 |877 |0.16 \n", + "|1980 |22 |814 |0.14 \n", + "|1981 |21 |362 |0.06 \n", + "|1982 |26 |764 |0.13 \n", + "|1983 |20 |809 |0.13 \n", + "|1984 |16 |223 |0.03 \n", + "|1985 |22 |1066 |0.15 \n", + "\n", + "+ Table 2.2 Worldwide airline fatalities, 1976–1985.\n", + "+ Death rate is passenger deaths per 100 million passenger miles.\n", + "+ Source: Statistical Abstract of the United States.\n", + "\n", + "Gelman, Andrew; Carlin, John B.; Stern, Hal S.; Dunson, David B.; Vehtari, Aki; Rubin, Donald B.. Bayesian Data Analysis, Third Edition (Chapman & Hall/CRC Texts in Statistical Science) (Page 60). CRC Press. Kindle Edition. \n" + ] + }, + { + "cell_type": "code", + "execution_count": 77, + "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", + "
DeathsFatalRateyearMiles
0734240.1919763863.0
1516250.1219774300.0
2754310.1519785027.0
3877310.1619795481.0
4814220.1419805814.0
5362210.0619816033.0
6764260.1319825877.0
7809200.1319836223.0
8223160.0319847433.0
91066220.1519857107.0
\n", + "
" + ], + "text/plain": [ + " Deaths Fatal Rate year Miles\n", + "0 734 24 0.19 1976 3863.0\n", + "1 516 25 0.12 1977 4300.0\n", + "2 754 31 0.15 1978 5027.0\n", + "3 877 31 0.16 1979 5481.0\n", + "4 814 22 0.14 1980 5814.0\n", + "5 362 21 0.06 1981 6033.0\n", + "6 764 26 0.13 1982 5877.0\n", + "7 809 20 0.13 1983 6223.0\n", + "8 223 16 0.03 1984 7433.0\n", + "9 1066 22 0.15 1985 7107.0" + ] + }, + "execution_count": 77, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from scipy.stats import poisson\n", + "from sklearn.linear_model import LinearRegression\n", + "import pystan\n", + "airline_df=pd.DataFrame(dict({'year':[x for x in range(1976,1986)],'Fatal':[24,25,31,31,22,21,26,20,16,22],'Deaths':[734,516,754,877,814,362,764,809,223,1066],'Rate':[.19,.12,.15,.16,.14,.06,.13,.13,.03,.15]}))\n", + "airline_df.set_index('year')\n", + "airline_df['Miles']=np.round(airline_df['Deaths']/airline_df['Rate'],0)\n", + "airline_df" + ] + }, + { + "cell_type": "code", + "execution_count": 66, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_62649727c9442f9bc7cbfce75826d859 NOW.\n" + ] + } + ], + "source": [ + "stan_code='''\n", + "data {\n", + " int deaths[10];\n", + "}\n", + "parameters {\n", + " real theta ; \n", + "}\n", + "model {\n", + "\n", + " // no prior here, what should we use?\n", + " deaths~poisson(theta);\n", + "}\n", + "\n", + "'''\n", + "sm_simple=pystan.StanModel(model_code=stan_code)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 67, + "metadata": {}, + "outputs": [], + "source": [ + "deaths=sm_simple.sampling(data=dict({'deaths':airline_df['Deaths']}))" + ] + }, + { + "cell_type": "code", + "execution_count": 68, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Inference for Stan model: anon_model_6ace3ad0d872dac6795ff2d2317d4760.\n", + "4 chains, each with iter=2000; warmup=1000; thin=1; \n", + "post-warmup draws per chain=1000, total post-warmup draws=4000.\n", + "\n", + " mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat\n", + "theta 691.94 0.21 8.16 676.54 686.33 691.83 697.4 708.53 1553 1.0\n", + "lp__ 3.8e4 0.02 0.69 3.8e4 3.8e4 3.8e4 3.8e4 3.8e4 1551 1.0\n", + "\n", + "Samples were drawn using NUTS at Sun Apr 15 18:46:11 2018.\n", + "For each parameter, n_eff is a crude measure of effective sample size,\n", + "and Rhat is the potential scale reduction factor on split chains (at \n", + "convergence, Rhat=1).\n" + ] + } + ], + "source": [ + "print(answer)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "From the stan output reported above, the 95% interval for the poisson rate is (676.5,708.5)." + ] + }, + { + "cell_type": "code", + "execution_count": 69, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_5173738091a79a17032e9cb1d2d99cd3 NOW.\n" + ] + } + ], + "source": [ + "stan_code='''\n", + "data {\n", + " int deaths[10];\n", + " vector[10] miles;\n", + "}\n", + "parameters {\n", + " real theta ; \n", + "}\n", + "model {\n", + "\n", + " // no prior here, what should we use?\n", + " deaths~poisson(miles*theta);\n", + "}\n", + "\n", + "'''\n", + "sm_weights=pystan.StanModel(model_code=stan_code)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 71, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/jet08013/anaconda3/lib/python3.6/site-packages/pystan/misc.py:399: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n", + " elif np.issubdtype(np.asarray(v).dtype, float):\n" + ] + } + ], + "source": [ + "deaths=sm_weights.sampling(data=dict({'deaths':airline_df['Deaths'],'miles':airline_df['Miles']}))" + ] + }, + { + "cell_type": "code", + "execution_count": 72, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Inference for Stan model: anon_model_5173738091a79a17032e9cb1d2d99cd3.\n", + "4 chains, each with iter=2000; warmup=1000; thin=1; \n", + "post-warmup draws per chain=1000, total post-warmup draws=4000.\n", + "\n", + " mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat\n", + "theta 0.12 4.4e-5 1.5e-3 0.12 0.12 0.12 0.12 0.12 1144 1.01\n", + "lp__ 3.8e4 0.02 0.71 3.8e4 3.8e4 3.8e4 3.8e4 3.8e4 1640 1.0\n", + "\n", + "Samples were drawn using NUTS at Mon Apr 16 06:52:19 2018.\n", + "For each parameter, n_eff is a crude measure of effective sample size,\n", + "and Rhat is the potential scale reduction factor on split chains (at \n", + "convergence, Rhat=1).\n" + ] + } + ], + "source": [ + "\n", + "print(deaths)" + ] + }, + { + "cell_type": "code", + "execution_count": 73, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHb1JREFUeJzt3Xl0lPd97/H3F4lFgEEsYtGCAQNiMXYEY8dL4niXwY5NbCBOcxviuOXeNDdN6lwa46Y3TdOeJJeem6Tn3pPUJ2nr9PQmGQHGDrFDHC/N0tjxCLHZWAZvzIxYxCJWAVp+94/nJxAgI6GR9Izm+bzO0ZmZ3zzS/OY58HxmnueZz5hzDhERiZ4BYU9ARETCoQAQEYkoBYCISEQpAEREIkoBICISUQoAEZGIUgCIiESUAkBEJKIUACIiEZUf9gQuZuzYsW7y5MlhT0NEpF+prq7e75wr6my5rA6AyZMnk0gkwp6GiEi/YmbvdWU57QISEYkoBYCISEQpAEREIkoBICISUQoAEZGIyuqzgEQkOtbVpFm1oZa6hkaKCwtYUVnOooqSsKeV0xQAIhK6dTVpVq7dSmNTCwDphkZWrt0KoBDoRdoFJCKhW7Wh9szGv01jUwurNtSGNKNoUACISOjqGhovaVx6hgJAREJXXFhwSePSMxQAIhK6FZXlFAzMO2esYGAeKyrLQ5pRNOggsIiEru1Ar84C6lsKABHJCosqSrTB72PaBSQiElEKABGRiFIAiIhElAJARCSiFAAiIhGlABARiSgFgIhIRHUaAGb2z2a2z8y2tRsbbWbPmdkOfznKj5uZ/aOZ7TSzLWY2r93vLPPL7zCzZb3zdEREpKu68g7gX4G7zht7FHjeOTcdeN7fBlgATPc/y4HvQRAYwFeBDwLXAl9tCw0REQlHpwHgnPs1cPC84fuAJ/z1J4BF7cZ/5AIvA4VmNhGoBJ5zzh10zh0CnuPCUBERkT7U3WMA451zuwH85Tg/XgIk2y2X8mPvNy4iIiHp6YPA1sGYu8j4hX/AbLmZJcwsUV9f36OTExGRs7obAHv9rh385T4/ngLK2i1XCtRdZPwCzrnHnXMx51ysqKiom9MTEZHOdDcAngbazuRZBjzVbvxT/myg64DDfhfRBuBOMxvlD/7e6cdERCQkndZBm9mPgZuBsWaWIjib55tA3MweBnYBS/zizwALgZ3ACeAhAOfcQTP7OvCqX+5vnXPnH1gWEZE+ZM51uCs+K8RiMZdIJMKehohIv2Jm1c65WGfL6ZPAIiIRpQAQEYkoBYCISETpO4FFRLLIm3uPUpVIUjZ6KJ+6fnKvPpYCQEQkZEdONvGzzXXEEyk2JxvIH2C9vvEHBYCISChaWx0vv3OAqkSKZ7ft5mRTKzPGD+crd8/iYxUljBk+uNfnoAAQEelD6YZG1lSnqKpOkjzYyGWD83lgXilLY2VcVToSs46ac3qHAkBEpJedbGrhudf3Ek8k+e3O/TgHN04bw5fuKKdyzgQKBuWFMi8FgIhIL3DO8VrdEeKJJE9tquNwYxMlhQX8+a3TWTy/lLLRQ8OeogJARKQnHTp+mnWb0sQTKbbvPsKg/AHcNWcCS2Nl3HDFGAYM6LtdPJ1RAIiIZKil1fGbHfVUJVI89/peTre0MrdkJF+/bw73Xl3CyKEDw55ihxQAIiLd9O7+46yuTrG6OsWeIycZNXQg/+W6y1kSK2XWxBFhT69TCgARkUtw4nQzz2zdQzyR5A/vHGSAwUdmFPHVj87m1lnjGJwfzgHd7lAAiIh0wjnHxl0NVCWSrN+ym2Onmpk8ZigrKst5YF4pE0YOCXuK3aIAEBF5H/uOnuTJjWniiSRv1R+nYGAed181kaWxMq6ZPKpPz9nvDQoAEZF2mlpaefGNfcQTKV6s3UdLqyN2+Si+9cBU7r6qmOGDc2ezmTvPREQkAzv2HiWeSPJkTZr9x05TdNlg/vTDU1kSK+WKouFhT69XKABEJLKOnGxi/ebdxBNJNvkStttmjWPJ/DJuLi8iPy+3G/MVACISKa2tjlfeOUhVIskz55WwLaooYWwflLBlCwWAiERC3ZkSthS7Dp4ItYQtWygARCRntZWwVVWn+M2OepyDG64YwyN3zAi1hC1bKACkX1lXk2bVhlrqGhopLixgRWU5iypKwp6WZJlt6cNUJZKs8yVsxSOH8Plbp7MkS0rYsoUCQPqNdTVpVq7dSmNTCxD0qq9cuxVAISAcOn6ap3wJ2+u+hK1yzgSWxkq54Yqx5GVRCVu2UABIv7FqQ+2ZjX+bxqYWVm2oVQBEVEur47c79xNPJHnutf5TwpYtFADSb9Q1NF7SuOSu9w6cLWHbfTgoYfvkdZNYMr+M2cXZX8KWLRQA0m8UFxaQ7mBjX1xYEMJspK+dON3Ms1v3UFWd5OW3gxK2m2YU8df3zOa2flbCli0UANJvrKgsP+cYAEDBwDxWVJaHOCvpTc45apJBCdvPNudOCVu2UABIv9G2n19nAeW++qOneLImRTyRYue+YzlXwpYtFADSryyqKNEGP0c1tbTyUm098USSF9/YR3OrY/7lo/jWA3NzroQtW2iNikiodu47SlUixZqNafYfO8XY4YN5+MNTWDK/jGnjcrOELVsoAESkzx092cT6LUEJW82usyVsS2NlfGRG7pewZQsFgIj0CeeCErZ4IskzW4MStunjolnCli0UACLSq3YfPlvC9t6BoITtfl/CdnVES9iyhQJARHrcqeYWfvX6PuKJJL/2JWzXTx3DF2+fzl1zJka+hC1bKABEpMe8VneYqkSKdZvSNJzwJWy3TGPx/DImjVEJW7ZRAIi0o7bRS9dw4jRPbaojnkjyWt0RBuUN4M4541kaK+PGaSphy2YKABFPbaNd11bCVpVI8ktfwnZlyQj+9r453Ht1MYVDB4U9RekCBYCIF2bbaH9557HrwAmqqpNnStgKhw7kjz44iSWxUuYUjwx7enKJMgoAM/sL4E8AB2wFHgImAj8BRgMbgT92zp02s8HAj4D5wAHg4865dzN5fJGeFFbbaLa/82g83cKz24Jz9l9++yBmcNP0Ir5y92xun60Stv6s2wFgZiXAnwOznXONZhYHHgQWAt92zv3EzL4PPAx8z18ecs5NM7MHgW8BH8/4GYj0kLDaRrPxew6cc2xKNhBPpFi/uY6jp5qZNHoo/+POGdw/r1QNrDki011A+UCBmTUBQ4HdwK3AH/n7nwD+hiAA7vPXAVYD/8fMzDnnMpyDSI8Iq200m77noP7oKdbVpIknkuzYd4whAwewcG5Qwnbt5NEM0AHdnNLtAHDOpc3sH4BdQCPwS6AaaHDONfvFUkDbS5gSIOl/t9nMDgNjgP3dnYNITwqrbTTs7znoqIStYlIh37h/LvdcNZHLhuhbtXJVJruARhG8qp8CNABVwIIOFm17hd/RS4cLXv2b2XJgOcCkSZO6Oz2RbgmjbTSsdx4dlbB95kNTWDK/lOnjL+vVx5bskMkuoNuBd5xz9QBmtha4ASg0s3z/LqAUqPPLp4AyIGVm+cBI4OD5f9Q59zjwOEAsFtPuIcl5ffnO4+jJJn7uS9g27mogb4Bx68yghO3m8iIGqoQtUjIJgF3AdWY2lGAX0G1AAngRWExwJtAy4Cm//NP+9u/9/S9o/79IoDffebSVsFUlUjyzdTeNTS1MGzecxxbOZFFFCeMu07dqRVUmxwBeMbPVBKd6NgM1BK/cfw78xMz+zo/90P/KD4F/M7OdBK/8H8xk4iJyceeXsA0fnM+iimKWxMqoKCtUCZtg2fwiPBaLuUQiEfY0RPqN9iVsv9lRT6uD66aOZmmsjLuunMDQQfrsZxSYWbVzLtbZcvrXINKJ/vAp3dfrjhBPJM+UsE0cOYTP3TKNxfNLuXzMsLCnJ1lKASByEdn8Kd22Eraq6iTb0kEJ2x2+hO1DKmGTLlAASJf0h1fBvSHbPqXb0ur43c79xNuVsM0pHsHX7g1K2EYNUwmbdJ0CQDqVza+Ce1u2fEp314ETrPYlbHWHTzKyIChhWzy/lCtLVMIm3aMAkE5l26vgvhTmp3QbT7fwi9d2E381xe/fPoAZfHh6EY/dPYvbZ41nyECVsElmFADSqWx5FRyGvv6UrnOOzanDxBNJfrbpbAnbl+6YwQPzVcImPUsBIJ0Ku6smTH31Kd39x07x5MY0VdVJ3tyrEjbpGzkZAFE9YNlbwuqqyRa99Snd5nYlbC/4ErYPlKmETfpOzgVAlA9Y9pawWjJz1c59x6iqTrJ2Y5r6o6cYO3yQStgkFDkXAFE+YNmbwmjJzCXHTjXz8y11xBMpqt87RN4A45bycSyNlXLLzHEqYZNQ5FwARPmApWQX5xyvvnuIeCLJz7cEJWxXFA1TCZtkjZwLgCgfsJTssOfwSdZsTFGVSPKuStgki+VcAET9gKWE41RzC89vD0rYfv1mUML2wSmj+fyt01kwVyVskp1y7l+lDlhKX9q+25ew1aQ5dKKJCSOG8Gc3ByVsk8eqhE2yW84FAOiApfSuwyeaeHpzmngixdb0YZWwSb+VkwEg0tNaWx2/e2s/8USKDa/t4XRzK7MmjuBvPjqb+z5QohI26ZcUACIXkTx4gqrqFGuqU6QbGhlZMJBPXFPGkliZStik31MAiJznZFMLv9i2h3giyX++FZSwfWjaWB5dMJM7ZquETXKHAkCE4Jz9Lb6E7enNdRw92UzZ6AIe8SVsJTqNWHKQAkAi7cCxUzxZk6YqkaJ271GGDBzAgisnsiRWynVTxqiETXKaAkAip7mllf94Myhhe3772RK2v//YlXz06mJGqIRNIkIBIJHxVv0xqhIp1mxMUX/0FGOGDeKhGyezJFbGDJWwSQQpACSnHTvVzDNbdhNPJEmcKWErYkmsjFtVwiYRpwCQnNNWwlaVSPLzrbs5cbqFqUXDWLlgJh+bpxI2kTYKAMkZbSVsq6tTvLP/OMMG5XHv1UEJ27xJKmETOZ8CQPq1082tPL99L/FEkv/wJWzXThnN526ZxkKVsIlclP53SL/0xp4jxF9NsW5TmoPHTzNhxBA+e/MVLJ5fxhSVsIl0iQJA+o22Eraq6hRbUocZmGfcOXsCS2KlfHh6kUrYRC6RAkCyWmur4z/fOkA8keQXvoRt5oTL+KovYRutEjaRblMASFZKHjzB6urggG66oZERQ/J58JoylsbKmFM8Qgd0RXqAAkCyxsmmFja8FpSw/W7n2RK2Ly+YyZ0qYRPpcQoACZVK2ETCowCQUJxfwjY4fwAL56qETaQvKQCkzzS3tPLrHfXEX03xq+17aW51XK0SNpHQKACk171df+zMt2rtUwmbSNZQAEivaCthq6pO8uq755aw3VI+jkH5KmETCZsCQHqMc47Ee4eIv3puCdujC2Zyf0UJ40aohE0kmygAJGN7jwQlbFWJsyVsH72qmKXXlDJv0iidsy+SpTIKADMrBH4AXAk44DNALfBTYDLwLrDUOXfIgq3Ad4GFwAng0865jZk8voTndHMrL7yxl3gixUu1+84pYVtw5QSGDdZrC5Fsl+n/0u8Cv3DOLTazQcBQ4DHgeefcN83sUeBR4MvAAmC6//kg8D1/Kf3IG3uOUJVI8WRNUMI2fsRglbCJ9FPdDgAzGwHcBHwawDl3GjhtZvcBN/vFngBeIgiA+4AfOecc8LKZFZrZROfc7m7PXvrE4cYmnt5cR1UieaaE7Y7Z41kSK+MmlbCJ9FuZvAOYCtQD/2JmVwPVwBeA8W0bdefcbjMb55cvAZLtfj/lxxQAWai11fH7t30J27Y9nPIlbP/zntksqlAJm0guyCQA8oF5wOedc6+Y2XcJdve8n45eJroLFjJbDiwHmDRpUgbTk+5IHQpK2KoSZ0vYPq4SNpGclEkApICUc+4Vf3s1QQDsbdu1Y2YTgX3tli9r9/ulQN35f9Q59zjwOEAsFrsgIKTntZWwVSVS/O6t/YBK2ESioNsB4JzbY2ZJMyt3ztUCtwGv+59lwDf95VP+V54G/ruZ/YTg4O9h7f8Pj3OObekjxBNJntqU5sjJZkpHFfDF22bwwPwSSkcNDXuKItLLMj0L6PPAv/szgN4GHgIGAHEzexjYBSzxyz5DcAroToLTQB/K8LGlGw4cO8W6TcEB3Tf2BCVsC66cwJJYGddPVQmbSJRkFADOuU1ArIO7butgWQd8LpPHk+5pbmnlNzv2E08k+dX2vTS1BCVsf7coKGEbWaASNpEo0qd1ctg7+49TlUiyZmOKvUeCErZl1wclbOUTVMImEnUKgBxz/FQzz2zdTVUixR/ePcgAg1vKx/G1e8u4daZK2ETkLAVADnDOUf3eIeKJJOu3+BK2scP48l0zuX9eCeNVwiYiHVAA9GP7jpxkzcY0VYkkb/sStnuumsjSWBnzL1cJm4hcnAKgnwlK2PZRlUjy0pv1tLQ6rp08ms/efAUL505UCZuIdJm2Fv1E7Z6jVCWSPFmT5sDx04y7bDD/9aapLJ5fytSi4WFPT0T6IQVAFjvc2MTPfAnbZl/Cdvus8SyNlfHh6WPJz9MBXRHpPgVAlmltdbzsS9iebVfC9tf3zGbRB4oZM3xw2FMUkRyhAMgSqUMnWFOdpqo6SepQI5cNyWdJrJSlsTLmlozUAV0R6XEKgBCdbGrhl6/vpSqR5Lc79+NcUMK2orKcyjkTVMImIr1KAdDHOiphKyks4Au3TeeBeaWUjVYJm4j0DQVAL1hXk2bVhlrqGhopLixgRWU5N80oYl1NmrgvYRvkS9iWqoRNREKiAOhh62rSrFy7lcamFgDSDY08Et+EmdHS6riqdCRfX3Ql915VzMihKmETkfAoAHrYqg21Zzb+bVodDBs0gDWfvYGZE0aENDMRkXMpAHpIWwlbuqGxw/tPnGrRxl9EsooCIAPOOTbuOkT81RTrt9Rx/HQLeQOCXT3nKy4sCGGGIiLvTwHQDfuOnGStP6D7dv1xhg7K4+65E1l6TRmpgyd47Mlt5+wGKhiYx4rK8hBnLCJyIQVAFzW1nC1he7E2KGG7ZvIo/ttHruDudiVs10wejZldcBbQooqSkJ+BiMi5FACdeHPvUeKvnlvCttyXsF3xPiVsiypKtMEXkaynAOjAkZNBCVs8kWJzsoH8Ab6E7ZpSbppepBI2EckJCgCvtdXx8jsHqEqkeHbbbk42tTJj/HC+cvcsPlZRohI2Eck5kQ+AdEMja6pTVFUnSR5s5LLB+TwwLyhhu6pUJWwikrsiHQBv7DnCgu/+Bufgxmlj+NIdQQlbwSCVsIlI7ot0AJSPv4zHFszirisnqIRNRCIn0gFgZvzpTVPDnoaISCh0OouISEQpAEREIkoBICISUQoAEZGIUgCIiESUAkBEJKIUACIiEaUAEBGJKAWAiEhEKQBERCJKASAiElEKABGRiFIAiIhEVMYBYGZ5ZlZjZuv97Slm9oqZ7TCzn5rZID8+2N/e6e+fnOlji4hI9/XEO4AvANvb3f4W8G3n3HTgEPCwH38YOOScmwZ82y8nIiIhySgAzKwUuBv4gb9twK3Aar/IE8Aif/0+fxt//22m71sUEQlNpu8AvgP8JdDqb48BGpxzzf52Cijx10uAJIC//7Bf/hxmttzMEmaWqK+vz3B6IiLyfrodAGZ2D7DPOVfdfriDRV0X7js74NzjzrmYcy5WVFTU3emJiEgnMvlKyBuBe81sITAEGEHwjqDQzPL9q/xSoM4vnwLKgJSZ5QMjgYMZPL6IiGSg2+8AnHMrnXOlzrnJwIPAC865TwIvAov9YsuAp/z1p/1t/P0vOOcueAcgIiJ9ozc+B/Bl4BEz20mwj/+HfvyHwBg//gjwaC88toiIdFEmu4DOcM69BLzkr78NXNvBMieBJT3xeCIikrkeCQARCayrSbNqQy11DY0UFxaworKcRRUlnf+iSAgUACI9ZF1NmpVrt9LY1AJAuqGRlWu3AigEJCupC0ikh6zaUHtm49+msamFVRtqQ5qRyMUpAER6SF1D4yWNi4RNASDSQ4oLCy5pXCRsCgCRHrKispyCgXnnjBUMzGNFZXlIMxK5OB0EFukhbQd6dRaQ9BcKAJEetKiiRBt86Te0C0hEJKIUACIiEaUAEBGJKAWAiEhEKQBERCJKASAiElEKABGRiFIAiIhElAJARCSiFAAiIhGlABARiSgFgIhIRCkAREQiSgEgIhJRCgARkYhSAIiIRJQCQEQkohQAIiIRpQAQEYkoBYCISEQpAEREIkoBICISUQoAEZGIUgCIiESUAkBEJKLyw56AiEhfWFeTZtWGWuoaGikuLGBFZTmLKkrCnlaoFAAikvPW1aRZuXYrjU0tAKQbGlm5ditApENAu4BEJOet2lB7ZuPfprGphVUbakOaUXZQAIhIzqtraLyk8ahQAIhIzisuLLik8ajodgCYWZmZvWhm283sNTP7gh8fbWbPmdkOfznKj5uZ/aOZ7TSzLWY2r6eehIjIxayoLKdgYN45YwUD81hRWR7SjLJDJu8AmoEvOedmAdcBnzOz2cCjwPPOuenA8/42wAJguv9ZDnwvg8cWEemyRRUlfOP+uZQUFmBASWEB37h/bqQPAEMGZwE553YDu/31o2a2HSgB7gNu9os9AbwEfNmP/8g554CXzazQzCb6vyMi0qsWVZREfoN/vh45BmBmk4EK4BVgfNtG3V+O84uVAMl2v5byY+f/reVmljCzRH19fU9MT0REOpBxAJjZcGAN8EXn3JGLLdrBmLtgwLnHnXMx51ysqKgo0+mJiMj7yCgAzGwgwcb/351za/3wXjOb6O+fCOzz4ymgrN2vlwJ1mTy+iIh0XyZnARnwQ2C7c+5/t7vraWCZv74MeKrd+Kf82UDXAYe1/19EJDyZVEHcCPwxsNXMNvmxx4BvAnEzexjYBSzx9z0DLAR2AieAhzJ4bBERyZAFJ+VkJzOrB97rg4caC+zvg8fpz7SOukbrqXNaR12TyXq63DnX6UHUrA6AvmJmCedcLOx5ZDOto67Reuqc1lHX9MV6UhWEiEhEKQBERCJKARB4POwJ9ANaR12j9dQ5raOu6fX1pGMAIiIRpXcAIiIRldMBYGZ5ZlZjZuv97Slm9oqvqv6pmQ3y44P97Z3+/snt/sZKP15rZpXhPJPeYWbvmtlWM9tkZgk/dsl13ma2zC+/w8yWvd/j9Ve+uHC1mb3h68+v13o6l5mV+39HbT9HzOyLWk/nMrO/8PX528zsx2Y2JNTtknMuZ3+AR4D/B6z3t+PAg/7694HP+ut/BnzfX38Q+Km/PhvYDAwGpgBvAXlhP68eXD/vAmPPG/tfwKP++qPAt/z1hcCzBJ1O1wGv+PHRwNv+cpS/Pirs59bD6+kJ4E/89UFAodbTRddXHrAHuFzr6Zz1UgK8AxT423Hg02Ful0JfKb24sksJvo/gVmC9/4e2H8j3918PbPDXNwDX++v5fjkDVgIr2/3NM8vlws/7BEAtMNFfnwjU+uv/BHzi/OWATwD/1G78nOX6+w8wwv+nNa2nLq+zO4HfaT1dsF7aGpFH++3MeqAyzO1SLu8C+g7wl0Crvz0GaHDONfvb7euoz1RV+/sP++W7VGHdjzngl2ZWbWbL/dil1nnn+jqaCtQD/+J3J/7AzIah9XQxDwI/9te1njznXBr4B4KKnN0E25lqQtwu5WQAmNk9wD7nXHX74Q4WdZ3c16UK637sRufcPIJva/ucmd10kWWjuo7ygXnA95xzFcBxzn7LXUeiup4A8Puv7wWqOlu0g7GcXk/++Md9BLttioFhBP/3ztdn26WcDACCorp7zexd4CcEu4G+AxSaWVsBXvs66jNV1f7+kcBBcrzC2jlX5y/3AU8C13Lpdd45vY4Inl/KOfeKv72aIBC0njq2ANjonNvrb2s9nXU78I5zrt451wSsBW4gxO1STgaAc26lc67UOTeZ4O3oC865TwIvAov9YudXVbedbbDYL+/8+IP+aPwUgu8z/kMfPY1eZWbDzOyytusE+223cel13huAO81slH+Fc6cfywnOuT1A0szavj38NuB1tJ7ezyc4u/sHtJ7a2wVcZ2ZDzcw4+28pvO1S2AdG+uDAy82cPQtoql9ROwneog7240P87Z3+/qntfv+vCI6y1wILwn4+PbhephKcSbAZeA34Kz8+huDg+Q5/OdqPG/B//brYCsTa/a3P+HW3E3go7OfWC+vqA0AC2AKsIzg7RevpwvU0FDgAjGw3pvV07jr6GvAGwYutfyM4kye07ZI+CSwiElE5uQtIREQ6pwAQEYkoBYCISEQpAEREIkoBICISUQoAEZGIUgCIiESUAkBEJKL+P42HafhjjkW6AAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig,ax=plt.subplots(1)\n", + "ax.scatter(airline_df['Miles'],airline_df['Deaths'] )\n", + "ax.plot(np.linspace(4000,8000,10),(.12)*np.linspace(4000,8000,10))\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 90, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'copy_X': True, 'fit_intercept': True, 'n_jobs': 1, 'normalize': False}" + ] + }, + "execution_count": 90, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lr=LinearRegression()\n", + "lr.fit(airline_df['Miles'].values.reshape(-1,1),airline_df['Deaths'].values.reshape(-1,1))\n", + "deaths_pred=lr.predict(np.linspace(4000,8000,10).reshape(-1,1))\n", + "ax.plot(np.linspace(4000,8000,10),deaths_pred)\n", + "plt.show()\n", + "lr.get_params()" + ] + }, + { + "cell_type": "code", + "execution_count": 84, + "metadata": {}, + "outputs": [ + { + "ename": "AttributeError", + "evalue": "'LinearRegression' object has no attribute 'slope'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mlr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mslope\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[0;31mAttributeError\u001b[0m: 'LinearRegression' object has no attribute 'slope'" + ] + } + ], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.6.5" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/BDA 2.11.13.md b/BDA 2.11.13.md new file mode 100644 index 0000000..f9126c8 --- /dev/null +++ b/BDA 2.11.13.md @@ -0,0 +1,32 @@ +### Problem 2.10.13 + +Discrete data: Table 2.2 gives the number of fatal accidents and deaths on scheduled airline flights per year over a ten-year period. +We use these data as a numerical example for fitting discrete data models. + +1. Assume that the numbers of fatal accidents in each year ar e independent with a Poisson(θ) distribution. Set a prior distribution forθand determine the posterior distribution based on the data from 1976 through 1985. Under this model, give a 95% predictive interval for the number of fatal accidents in 1986. You can use the normal approximation to the gamma and Poisson or compute using simulation. +2. Assume that the numbers of fatal accidents in each year fo llow independent Poisson distributions with a constant rate and an exposure in each year proportional to the number of passenger miles flown. Set a prior distribution forθand determine the posterior distribution based on the data for 1976–1985. (Estimate the number of passenger miles flown in each year by dividing the appropriate columns of Table 2.2 and ignoring round-off errors.) Give a 95% predictive interval for the number of fatal iaccidents in 1986 under the assumption that 8 × 10 11 passenger miles are flown that year. +3. Repeat (1) above, replacing ‘fatal accidents’ with ‘passenger deaths.’ +4. Repeat (2) above, replacing ‘fatal accidents’ with ‘passenger deaths.’ +5. In which of the cases (1)–(4) above does the Poisson model seem more or less reasonable? Why? Discuss based on general principles, +without specific reference to the numbers in Table 2.2. Incidentally, in 1986, there were 22 fatal accidents, +546 passenger deaths, and a death rate of 0.06 per 100 million miles flown. We return to this example in Exercises 3.12, 6.2, 6.3, and 8.14. + +|Year |Fatal accidents |Passenger deaths |Death rate +|---|---|---|---| +|1976 | 24 | 734 | 0.19 +|1977 |25 |516 |0.12 +|1978 |31 |754 |0.15 +|1979 |31 |877 |0.16 +|1980 |22 |814 |0.14 +|1981 |21 |362 |0.06 +|1982 |26 |764 |0.13 +|1983 |20 |809 |0.13 +|1984 |16 |223 |0.03 +|1985 |22 |1066 |0.15 + +Table 2.2 Worldwide airline fatalities, 1976–1985. + +*Death rate is passenger deaths per 100 million passenger miles. +Source: Statistical Abstract of the United States.* + +Gelman, Andrew; Carlin, John B.; Stern, Hal S.; Dunson, David B.; Vehtari, Aki; Rubin, Donald B.. Bayesian Data Analysis, Third Edition (Chapman & Hall/CRC Texts in Statistical Science) (Page 60). CRC Press. Kindle Edition. diff --git a/BDA 2.11.13.txt b/BDA 2.11.13.txt deleted file mode 100644 index 37d948a..0000000 --- a/BDA 2.11.13.txt +++ /dev/null @@ -1,9 +0,0 @@ -13. Discrete data: Table 2.2 gives the number of fatal accidents and deaths on scheduled airline flights per year over a ten-year period. We use these data as a numerical example for fitting discrete data models. (a) Assume that the numbers of fatal accidents in each year ar e independent with a Poisson(θ) distribution. Set a prior distribution forθand determine the posterior distribution based on the data from 1976 through 1985. Under this model, give a 95% predictive interval for the number of fatal accidents in 1986. You can use the normal approximation to the gamma and Poisson or compute using simulation. (b) Assume that the numbers of fatal accidents in each year fo llow independent Poisson distributions with a constant rate and an exposure in each year proportional to the number of passenger miles flown. Set a prior distribution forθand determine the posterior distribution based on the data for 1976–1985. (Estimate the number of passenger miles flown in each year by dividing the appropriate columns of Table 2.2 and ignoring round-off errors.) Give a 95% predictive interval for the number of fatal iaccidents in 1986 under the assumption that 8 × 10 11 passenger miles are flown that year. (c) Repeat (a) above, replacing ‘fatal accidents’ with ‘passenger deaths.’ (d) Repeat (b) above, replacing ‘fatal accidents’ with ‘passenger deaths.’ (e) In which of the cases (a)–(d) above does the Poisson model seem more or less reasonable? Why? Discuss based on general principles, without specific reference to the numbers in Table 2.2. Incidentally, in 1986, there were 22 fatal accidents, 546 passenger deaths, and a death rate of 0.06 per 100 million miles flown. We return to this example in Exercises 3.12, 6.2, 6.3, and 8.14. - -Year Fatal Passenger Death accidents deaths rate 1976 24 734 0.19 1977 25 516 0.12 1978 31 754 0.15 1979 31 877 0.16 1980 22 814 0.14 1981 21 362 0.06 1982 26 764 0.13 1983 20 809 0.13 1984 16 223 0.03 1985 22 1066 0.15 Table 2.2 Worldwide airline fatalities, 1976–1985. Death rate is passen ger deaths per 100 million passenger miles. Source: Statistical Abstract of the United States. - -Gelman, Andrew; Carlin, John B.; Stern, Hal S.; Dunson, David B.; Vehtari, Aki; Rubin, Donald B.. Bayesian Data Analysis, Third Edition (Chapman & Hall/CRC Texts in Statistical Science) (Page 60). CRC Press. Kindle Edition. - -Gelman, Andrew; Carlin, John B.; Stern, Hal S.; Dunson, David B.; Vehtari, Aki; Rubin, Donald B.. Bayesian Data Analysis, Third Edition (Chapman & Hall/CRC Texts in Statistical Science) (Page 60). CRC Press. Kindle Edition. - -Gelman, Andrew; Carlin, John B.; Stern, Hal S.; Dunson, David B.; Vehtari, Aki; Rubin, Donald B.. Bayesian Data Analysis, Third Edition (Chapman & Hall/CRC Texts in Statistical Science) (Page 59). CRC Press. Kindle Edition.