week08b-poisson-regression-tropical-storms.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Poisson regression example\n",
    "\n",
    "## Modeling tropical storm counts\n",
    "\n",
    "The number of tropical storms is an example of a variable that is not expected to be normally distributed. Counts cannot be negative, are often skewed and tend to follow a Poisson distribution.\n",
    "\n",
    "This example follows the general approach of [Villarini et al. (2010)](https://journals.ametsoc.org/view/journals/mwre/138/7/2010mwr3315.1.xml) in using Poisson regression to model the number of tropical storms based on climate indices.\n",
    "\n",
    "The goal of Poisson regression is to model a rate of occurrence, $\\Lambda_i$, where the rate changes for each observation $i$. This model can be expressed as a generalized linear model (GLM) with a *logarithmic link function*.\n",
    "\n",
    "$$\\log(\\Lambda_i) = \\beta_0 + \\beta_1 x_{1i} + \\beta_2 x_{2i} + ... + \\beta_n x_{ki}$$\n",
    "\n",
    "which is the same as\n",
    "\n",
    "$$\\Lambda_i = \\exp(\\beta_0 + \\beta_1 x_{1i} + \\beta_2 x_{2i} + ... + \\beta_n x_{ki})$$\n",
    "\n",
    "This is a model for the rate of occurrence $\\Lambda_i$ as a function of $k$ predictor variables. In this example, the rate of occurrence $\\Lambda_i$ is the number of storm counts per year and each climate index is a predictor variable. \n",
    "\n",
    "* The logarithmic link function is useful because $\\log(\\Lambda_i)$ can be positive or negative, but $\\Lambda_i$ must be positive\n",
    "\n",
    "* If the data being modeled is a standard Poisson random variable, the model simplifies to $\\Lambda_i = \\exp({\\beta_0})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy import stats\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use a special function to load the climate indices from NOAA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_psl_file(psl_file):\n",
    "    '''\n",
    "    Read a data file in in NOAA Physical Sciences Laboratory (PSL) format\n",
    "    \n",
    "    Input: String containing the path to a PSL data file\n",
    "    Returns: Pandas dataframe with monthly data in columns and the year as the index\n",
    "    \n",
    "    For a description of the PSL format see: https://psl.noaa.gov/gcos_wgsp/Timeseries/standard/\n",
    "    '''\n",
    "    \n",
    "    f = open(psl_file, \"r\")\n",
    "    all_lines = f.readlines()\n",
    "    start_year = all_lines[0].split()[0]\n",
    "    end_year = all_lines[0].split()[1]\n",
    "\n",
    "    for i in range(1,len(all_lines)):\n",
    "        stri = all_lines[i].find(end_year)\n",
    "        if stri>=0:\n",
    "            end_index = i\n",
    "\n",
    "    missing_val = float(all_lines[end_index+1])\n",
    "    nrows = int(end_year)-int(start_year)+1\n",
    "    df = pd.read_csv(psl_file,skiprows=1,nrows=nrows,sep='\\s+',header=None,na_values=missing_val)\n",
    "    df = df.rename(columns={0:'year'})\n",
    "    df = df.set_index('year',drop=True)\n",
    "    \n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfsoi = read_psl_file('data/tropical-storms/soi.data')\n",
    "dftna = read_psl_file('data/tropical-storms/tna.data')\n",
    "dfnao = read_psl_file('data/tropical-storms/nao.data')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the tropical storm data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "dftrop = pd.read_csv('data/tropical-storms/tropical.txt',sep='\\t')\n",
    "dftrop = dftrop.set_index('Year',drop=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tropical storms do not happen in winter, so we average the climate indices during the relevant months."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "dftrop['SOI'] = dfsoi.loc[:,5:6].mean(axis=1)\n",
    "dftrop['TNA'] = dftna.loc[:,5:6].mean(axis=1)\n",
    "dftrop['NAO'] = dfnao.loc[:,5:6].mean(axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercises\n",
    "* How would you find the NAO averaged over July-September?\n",
    "* How would you find the NAO for years 2010-2020, averaged over July-September?\n",
    "* How would you use `dfnao.iloc` to achieve the same result?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drop all rows with any NaN values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = dftrop.dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Year</th>\n",
       "      <th>NamedStorms</th>\n",
       "      <th>NamedStormDays</th>\n",
       "      <th>Hurricanes</th>\n",
       "      <th>HurricaneDays</th>\n",
       "      <th>MajorHurricanes</th>\n",
       "      <th>MajorHurricaneDays</th>\n",
       "      <th>AccumulatedCycloneEnergy</th>\n",
       "      <th>SOI</th>\n",
       "      <th>TNA</th>\n",
       "      <th>NAO</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Year</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1951</th>\n",
       "      <td>1951</td>\n",
       "      <td>12</td>\n",
       "      <td>67.00</td>\n",
       "      <td>8</td>\n",
       "      <td>34.25</td>\n",
       "      <td>3</td>\n",
       "      <td>4.50</td>\n",
       "      <td>126.3</td>\n",
       "      <td>-0.40</td>\n",
       "      <td>0.155</td>\n",
       "      <td>-0.925</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1952</th>\n",
       "      <td>1952</td>\n",
       "      <td>11</td>\n",
       "      <td>45.75</td>\n",
       "      <td>5</td>\n",
       "      <td>15.50</td>\n",
       "      <td>2</td>\n",
       "      <td>2.50</td>\n",
       "      <td>69.1</td>\n",
       "      <td>1.20</td>\n",
       "      <td>0.205</td>\n",
       "      <td>-0.545</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1953</th>\n",
       "      <td>1953</td>\n",
       "      <td>14</td>\n",
       "      <td>61.75</td>\n",
       "      <td>7</td>\n",
       "      <td>19.00</td>\n",
       "      <td>3</td>\n",
       "      <td>5.00</td>\n",
       "      <td>98.5</td>\n",
       "      <td>-1.30</td>\n",
       "      <td>0.120</td>\n",
       "      <td>0.425</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1954</th>\n",
       "      <td>1954</td>\n",
       "      <td>16</td>\n",
       "      <td>57.00</td>\n",
       "      <td>7</td>\n",
       "      <td>26.00</td>\n",
       "      <td>3</td>\n",
       "      <td>7.00</td>\n",
       "      <td>104.4</td>\n",
       "      <td>0.50</td>\n",
       "      <td>0.045</td>\n",
       "      <td>0.015</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1955</th>\n",
       "      <td>1955</td>\n",
       "      <td>13</td>\n",
       "      <td>85.50</td>\n",
       "      <td>9</td>\n",
       "      <td>43.00</td>\n",
       "      <td>4</td>\n",
       "      <td>8.50</td>\n",
       "      <td>164.7</td>\n",
       "      <td>1.95</td>\n",
       "      <td>-0.065</td>\n",
       "      <td>-0.530</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2016</th>\n",
       "      <td>2016</td>\n",
       "      <td>15</td>\n",
       "      <td>81.00</td>\n",
       "      <td>7</td>\n",
       "      <td>27.75</td>\n",
       "      <td>4</td>\n",
       "      <td>10.25</td>\n",
       "      <td>141.3</td>\n",
       "      <td>0.90</td>\n",
       "      <td>0.385</td>\n",
       "      <td>-0.400</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017</th>\n",
       "      <td>2017</td>\n",
       "      <td>17</td>\n",
       "      <td>93.00</td>\n",
       "      <td>10</td>\n",
       "      <td>51.75</td>\n",
       "      <td>6</td>\n",
       "      <td>19.25</td>\n",
       "      <td>224.9</td>\n",
       "      <td>-0.15</td>\n",
       "      <td>0.590</td>\n",
       "      <td>-0.685</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2018</th>\n",
       "      <td>2018</td>\n",
       "      <td>15</td>\n",
       "      <td>87.25</td>\n",
       "      <td>8</td>\n",
       "      <td>26.75</td>\n",
       "      <td>2</td>\n",
       "      <td>5.00</td>\n",
       "      <td>129.0</td>\n",
       "      <td>0.20</td>\n",
       "      <td>-0.440</td>\n",
       "      <td>1.715</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019</th>\n",
       "      <td>2019</td>\n",
       "      <td>18</td>\n",
       "      <td>68.50</td>\n",
       "      <td>6</td>\n",
       "      <td>23.25</td>\n",
       "      <td>3</td>\n",
       "      <td>10.00</td>\n",
       "      <td>129.8</td>\n",
       "      <td>-0.70</td>\n",
       "      <td>0.210</td>\n",
       "      <td>-1.585</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020</th>\n",
       "      <td>2020</td>\n",
       "      <td>30</td>\n",
       "      <td>118.00</td>\n",
       "      <td>13</td>\n",
       "      <td>34.75</td>\n",
       "      <td>6</td>\n",
       "      <td>8.75</td>\n",
       "      <td>179.8</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.620</td>\n",
       "      <td>-0.085</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>70 rows × 11 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      Year  NamedStorms  NamedStormDays  Hurricanes  HurricaneDays  \\\n",
       "Year                                                                 \n",
       "1951  1951           12           67.00           8          34.25   \n",
       "1952  1952           11           45.75           5          15.50   \n",
       "1953  1953           14           61.75           7          19.00   \n",
       "1954  1954           16           57.00           7          26.00   \n",
       "1955  1955           13           85.50           9          43.00   \n",
       "...    ...          ...             ...         ...            ...   \n",
       "2016  2016           15           81.00           7          27.75   \n",
       "2017  2017           17           93.00          10          51.75   \n",
       "2018  2018           15           87.25           8          26.75   \n",
       "2019  2019           18           68.50           6          23.25   \n",
       "2020  2020           30          118.00          13          34.75   \n",
       "\n",
       "      MajorHurricanes  MajorHurricaneDays  AccumulatedCycloneEnergy   SOI  \\\n",
       "Year                                                                        \n",
       "1951                3                4.50                     126.3 -0.40   \n",
       "1952                2                2.50                      69.1  1.20   \n",
       "1953                3                5.00                      98.5 -1.30   \n",
       "1954                3                7.00                     104.4  0.50   \n",
       "1955                4                8.50                     164.7  1.95   \n",
       "...               ...                 ...                       ...   ...   \n",
       "2016                4               10.25                     141.3  0.90   \n",
       "2017                6               19.25                     224.9 -0.15   \n",
       "2018                2                5.00                     129.0  0.20   \n",
       "2019                3               10.00                     129.8 -0.70   \n",
       "2020                6                8.75                     179.8  0.05   \n",
       "\n",
       "        TNA    NAO  \n",
       "Year                \n",
       "1951  0.155 -0.925  \n",
       "1952  0.205 -0.545  \n",
       "1953  0.120  0.425  \n",
       "1954  0.045  0.015  \n",
       "1955 -0.065 -0.530  \n",
       "...     ...    ...  \n",
       "2016  0.385 -0.400  \n",
       "2017  0.590 -0.685  \n",
       "2018 -0.440  1.715  \n",
       "2019  0.210 -1.585  \n",
       "2020  0.620 -0.085  \n",
       "\n",
       "[70 rows x 11 columns]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histogram of the named storm counts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAATAklEQVR4nO3df4xc13ne8e9TyUZUriI5UbyRKKV0A0GAIsaKtaBTuAh269igaKNKDDcVIThS4oB2ELUOoAJmUyBxWwRQ2spFAhlWmJqQDTTaGHDUCBJjWzCylQXYsUmBNqlKqlmFaUkaJBTZtNcWYNB6+8deotvVDHd3fnCXZ74fYDEz5557z3lxyWdnz8zcSVUhSWrX39voCUiSxsugl6TGGfSS1DiDXpIaZ9BLUuMu3+gJ9HLNNdfUtm3bxj7O9773PbZs2TL2cTYDa23PpNQJk1PrMHUeOnTopar6iV7bNmXQb9u2jYMHD459nIWFBWZnZ8c+zmZgre2ZlDphcmodps4kf9tvm0s3ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUuE35yVhdGrbtfWLDxj5+/7s2bGzpUuMzeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1btWLmiXZD7wbOFNVt3Rtfwbc1HW5Gvh2Vd3aY9/jwHeBHwLnqmpmJLOWJK3ZWq5e+TDwIPCp8w1V9c/P30/yAHD2AvvPVdVLg05QkjScVYO+qp5Ksq3XtiQBfgX4JyOelyRpRFJVq3daCvrHzy/dLGv/BeCj/ZZkkvwN8C2ggD+uqn0XGGMPsAdgenr6tvn5+bXWMLDFxUWmpqbGPs5mMI5aj5y80B9y47V961V9t03KeZ2UOmFyah2mzrm5uUP9snjYLx7ZDTxyge1vq6pTSd4IPJnk+ap6qlfH7pfAPoCZmZmanZ0dcmqrW1hY4GKMsxmMo9Z7NvKLR+6a7bttUs7rpNQJk1PruOoc+F03SS4H3gP8Wb8+VXWquz0DPArsGHQ8SdJghnl75S8Cz1fViV4bk2xJcuX5+8A7gaNDjCdJGsCqQZ/kEeBLwE1JTiR5f7fpTlYs2yS5LsmB7uE08HSSrwFfAZ6oqs+ObuqSpLVYy7tudvdpv6dH2ylgV3f/ReDNQ85PkjSkYV+MVWfbRr4wef+7NmxsSZufl0CQpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxq3lO2P3JzmT5Oiyto8kOZnkcPezq8++O5O8kORYkr2jnLgkaW3W8oz+YWBnj/b/XFW3dj8HVm5MchnwMeB24GZgd5Kbh5msJGn9Vg36qnoKeHmAY+8AjlXVi1X1A2AeuGOA40iShpCqWr1Tsg14vKpu6R5/BLgH+A5wELivqr61Yp/3Ajur6je6x+8D3lpV9/YZYw+wB2B6evq2+fn5wSpah8XFRaampkZyrCMnz47kOIPYvvWqVfuMstbzNmvN46h1M5qUOmFyah2mzrm5uUNVNdNr2+UDzufjwL8Hqrt9APj1FX3SY7++v1Wqah+wD2BmZqZmZ2cHnNraLSwsMKpx7tn7xEiOM4jjd82u2meUtZ63WWseR62b0aTUCZNT67jqHOhdN1V1uqp+WFWvAn/C0jLNSieAG5Y9vh44Nch4kqTBDRT0Sa5d9vCXgaM9un0VuDHJm5K8HrgTeGyQ8SRJg1t16SbJI8AscE2SE8DvAbNJbmVpKeY48IGu73XAf6mqXVV1Lsm9wOeAy4D9VfXsOIqQJPW3atBX1e4ezZ/o0/cUsGvZ4wPAa956KUm6ePxkrCQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxq0a9En2JzmT5Oiytv+Y5PkkX0/yaJKr++x7PMmRJIeTHBzhvCVJa7SWZ/QPAztXtD0J3FJVPwv8T+BfX2D/uaq6tapmBpuiJGkYqwZ9VT0FvLyi7fNVda57+GXg+jHMTZI0AqNYo/914C/7bCvg80kOJdkzgrEkSeuUqlq9U7INeLyqblnR/m+AGeA91eNASa6rqlNJ3sjScs+/6P5C6DXGHmAPwPT09G3z8/PrrWXdFhcXmZqaGsmxjpw8O5LjDGL71qtW7TPKWs/brDWPo9bNaFLqhMmpdZg65+bmDvVbIr980AkluRt4N/D2XiEPUFWnutszSR4FdgA9g76q9gH7AGZmZmp2dnbQqa3ZwsICoxrnnr1PjOQ4gzh+1+yqfUZZ63mbteZx1LoZTUqdMDm1jqvOgZZukuwEPgz806r6fp8+W5Jcef4+8E7gaK++kqTxWcvbKx8BvgTclOREkvcDDwJXAk92b518qOt7XZID3a7TwNNJvgZ8BXiiqj47liokSX2tunRTVbt7NH+iT99TwK7u/ovAm4eanSRpaAOv0Wvz2LaGtfL7tp/b0DV1SRvHSyBIUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWrcWr4cfH+SM0mOLmv7sSRPJvlGd/uGPvvuTPJCkmNJ9o5y4pKktVnLM/qHgZ0r2vYCX6iqG4EvdI//P0kuAz4G3A7cDOxOcvNQs5UkrduqQV9VTwEvr2i+A/hkd/+TwC/12HUHcKyqXqyqHwDz3X6SpIsoVbV6p2Qb8HhV3dI9/nZVXb1s+7eq6g0r9nkvsLOqfqN7/D7grVV1b58x9gB7AKanp2+bn58fqKD1WFxcZGpqaiTHOnLy7EiOMy7TV8DpVzZ6FqOzfetVfbeN8rxuZpNSJ0xOrcPUOTc3d6iqZnptu3yoWV1YerT1/a1SVfuAfQAzMzM1Ozs7pmn9PwsLC4xqnHv2PjGS44zLfdvP8cCRcZ7ui+v4XbN9t43yvG5mk1InTE6t46pz0HfdnE5yLUB3e6ZHnxPADcseXw+cGnA8SdKABg36x4C7u/t3A3/Ro89XgRuTvCnJ64E7u/0kSRfRWt5e+QjwJeCmJCeSvB+4H3hHkm8A7+gek+S6JAcAquoccC/wOeA54NNV9ex4ypAk9bPqom1V7e6z6e09+p4Cdi17fAA4MPDsJElD85OxktQ4g16SGmfQS1LjDHpJapxBL0mNa+ejkpoo2y7wSeT7tp8b2yeVj9//rrEcVxonn9FLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEDB32Sm5IcXvbznSS/vaLPbJKzy/r87tAzliSty8BXr6yqF4BbAZJcBpwEHu3R9YtV9e5Bx5EkDWdUSzdvB/5XVf3tiI4nSRqRVNXwB0n2A89U1YMr2meBzwAngFPAv6qqZ/scYw+wB2B6evq2+fn5oee1msXFRaampkZyrCMnz47kOOMyfQWcfmWjZ3FxjLPW7VuvGs+BBzDKf7+b3aTUOkydc3Nzh6pqpte2oYM+yetZCvGfqarTK7b9KPBqVS0m2QX8YVXduNoxZ2Zm6uDBg0PNay0WFhaYnZ0dybEu9EUYm8F928/xwJHJ+J6Zcda6mb54ZJT/fje7Sal1mDqT9A36USzd3M7Ss/nTKzdU1XeqarG7fwB4XZJrRjCmJGmNRhH0u4FHem1I8pNJ0t3f0Y33dyMYU5K0RkP9fZvk7wPvAD6wrO2DAFX1EPBe4DeTnANeAe6sUbwoIElas6GCvqq+D/z4iraHlt1/EHhw5X6SpItnMl6dk0Zko15030wvAuvS4yUQJKlxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXFDBX2S40mOJDmc5GCP7UnyR0mOJfl6krcMM54kaf1G8VWCc1X1Up9ttwM3dj9vBT7e3UqSLpJxL93cAXyqlnwZuDrJtWMeU5K0TKpq8J2TvwG+BRTwx1W1b8X2x4H7q+rp7vEXgA9XVa9lnj3AHoDp6enb5ufnB57XWi0uLjI1NTWSYx05eXYkxxmX6Svg9CsbPYuLo8Vat2+96jVto/z3u9lNSq3D1Dk3N3eoqmZ6bRt26eZtVXUqyRuBJ5M8X1VPLdueHvv0/M3S/ZLYBzAzM1Ozs7NDTm11CwsLjGqce/Y+MZLjjMt928/xwJFRrNRtfi3Wevyu2de0jfLf72Y3KbWOq86hlm6q6lR3ewZ4FNixossJ4IZlj68HTg0zpiRpfQYO+iRbklx5/j7wTuDoim6PAb/avfvm54GzVfXNgWcrSVq3Yf6+nQYeTXL+OH9aVZ9N8kGAqnoIOADsAo4B3wd+bbjpSpLWa+Cgr6oXgTf3aH9o2f0CfmvQMSRJw/OTsZLUOINekhpn0EtS4wx6SWqcQS9JjWvr44PAtnV8QvW+7ec2/SdaJWlYPqOXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuOG+XLwG5L8VZLnkjyb5EM9+swmOZvkcPfzu8NNV5K0XsNcvfIccF9VPZPkSuBQkier6n+s6PfFqnr3EONIkoYw8DP6qvpmVT3T3f8u8BywdVQTkySNxkjW6JNsA34O+Osem/9Rkq8l+cskPzOK8SRJa5eqGu4AyRTw34Hfr6o/X7HtR4FXq2oxyS7gD6vqxj7H2QPsAZienr5tfn5+oPkcOXl2zX2nr4DTrww0zCXHWi9t27de9Zq2xcVFpqamNmA2F9+k1DpMnXNzc4eqaqbXtqGCPsnrgMeBz1XVR9fQ/zgwU1UvXajfzMxMHTx4cKA5rfcbph440tyXbPVkrZe24/e/6zVtCwsLzM7OXvzJbIBJqXWYOpP0Dfph3nUT4BPAc/1CPslPdv1IsqMb7+8GHVOStH7DPO15G/A+4EiSw13b7wA/BVBVDwHvBX4zyTngFeDOGnatSJK0LgMHfVU9DWSVPg8CDw46hiRpeG0tZErSCKzntb5RenjnlrEc10sgSFLjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklq3FBBn2RnkheSHEuyt8f2JPmjbvvXk7xlmPEkSes3cNAnuQz4GHA7cDOwO8nNK7rdDtzY/ewBPj7oeJKkwQzzjH4HcKyqXqyqHwDzwB0r+twBfKqWfBm4Osm1Q4wpSVqny4fYdyvwf5Y9PgG8dQ19tgLfXHmwJHtYetYPsJjkhSHmtib/Eq4BXhr3OJuBtV7a8gc9m5ur8wImota5Pxiqzn/Qb8MwQZ8ebTVAn6XGqn3AviHms25JDlbVzMUcc6NYa3smpU6YnFrHVecwSzcngBuWPb4eODVAH0nSGA0T9F8FbkzypiSvB+4EHlvR5zHgV7t33/w8cLaqXrNsI0kan4GXbqrqXJJ7gc8BlwH7q+rZJB/stj8EHAB2AceA7wO/NvyUR+qiLhVtMGttz6TUCZNT61jqTFXPJXNJUiP8ZKwkNc6gl6TGTWzQJzme5EiSw0kObvR8RinJ/iRnkhxd1vZjSZ5M8o3u9g0bOcdR6FPnR5Kc7M7r4SS7NnKOo5DkhiR/leS5JM8m+VDX3uI57Vdri+f1R5J8JcnXulr/bdc+8vM6sWv0SY4DM1XV3IcwkvwCsMjSp5Jv6dr+A/ByVd3fXZfoDVX14Y2c57D61PkRYLGq/tNGzm2Uuk+TX1tVzyS5EjgE/BJwD+2d0361/grtndcAW6pqMcnrgKeBDwHvYcTndWKf0besqp4CXl7RfAfwye7+J1n6z3NJ61Nnc6rqm1X1THf/u8BzLH3CvMVz2q/W5nSXhlnsHr6u+ynGcF4nOegL+HySQ93lF1o3ff4zDN3tGzd4PuN0b3e11P0tLGcsl2Qb8HPAX9P4OV1RKzR4XpNcluQwcAZ4sqrGcl4nOejfVlVvYekKm7/VLQPo0vdx4KeBW1m6ptIDGzqbEUoyBXwG+O2q+s5Gz2ecetTa5Hmtqh9W1a0sXTVgR5JbxjHOxAZ9VZ3qbs8Aj7J0Nc6WnT5/5dDu9swGz2csqup095/nVeBPaOS8dmu4nwH+a1X9edfc5DntVWur5/W8qvo2sADsZAzndSKDPsmW7oUekmwB3gkcvfBel7zHgLu7+3cDf7GBcxmbFZfB/mUaOK/di3afAJ6rqo8u29TcOe1Xa6Pn9SeSXN3dvwL4ReB5xnBeJ/JdN0n+IUvP4mHpMhB/WlW/v4FTGqkkjwCzLF3a9TTwe8B/Az4N/BTwv4F/VlWX9AuZfeqcZenP+wKOAx+41K+vlOQfA18EjgCvds2/w9LadWvntF+tu2nvvP4sSy+2XsbSk+5PV9W/S/LjjPi8TmTQS9IkmcilG0maJAa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJatz/BY/6a5bqprMvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "df['NamedStorms'].hist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Poisson distribution - model as a constant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Generalized Linear Model Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>      <td>NamedStorms</td>   <th>  No. Observations:  </th>  <td>    70</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                  <td>GLM</td>       <th>  Df Residuals:      </th>  <td>    69</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model Family:</th>         <td>Poisson</td>     <th>  Df Model:          </th>  <td>     0</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Link Function:</th>          <td>log</td>       <th>  Scale:             </th> <td>  1.0000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -206.92</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Tue, 08 Mar 2022</td> <th>  Deviance:          </th> <td>  113.91</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>14:16:31</td>     <th>  Pearson chi2:      </th>  <td>  123.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>          <td>4</td>        <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>     <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>    2.4979</td> <td>    0.034</td> <td>   72.869</td> <td> 0.000</td> <td>    2.431</td> <td>    2.565</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                 Generalized Linear Model Regression Results                  \n",
       "==============================================================================\n",
       "Dep. Variable:            NamedStorms   No. Observations:                   70\n",
       "Model:                            GLM   Df Residuals:                       69\n",
       "Model Family:                 Poisson   Df Model:                            0\n",
       "Link Function:                    log   Scale:                          1.0000\n",
       "Method:                          IRLS   Log-Likelihood:                -206.92\n",
       "Date:                Tue, 08 Mar 2022   Deviance:                       113.91\n",
       "Time:                        14:16:31   Pearson chi2:                     123.\n",
       "No. Iterations:                     4                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept      2.4979      0.034     72.869      0.000       2.431       2.565\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = smf.glm(formula='NamedStorms ~ 1',\n",
    "                 data=df,\n",
    "                 family=sm.families.Poisson()).fit()\n",
    "result.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fb10199b4f0>]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAon0lEQVR4nO3deXxU5b3H8c/vTBIg7AqE3YAEAgIKsgbEqIAsKm21rVirtLdSqrZa8bbc2/a2t9fea73ValurRb1trQulVlsQZFGMG4ssohDCEhYlgCyCQNiSmfPcP86gISTkJJnJM3Pm93695pVZzpn5PpzMj5PnPOc5YoxBKaVUcDm2AyillIovLfRKKRVwWuiVUirgtNArpVTAaaFXSqmAS7MdoCpt2rQx2dnZdVr32LFjNG3aNLaBLAlKW4LSDtC2JKKgtAPq15bVq1cfMMa0req1hCz02dnZrFq1qk7rFhQUkJ+fH9tAlgSlLUFpB2hbElFQ2gH1a4uIfFjda9p1o5RSAaeFXimlAk4LvVJKBZwWeqWUCjgt9EopFXBa6JVSKuC00CulVMBpoVfJR6fWVqpWEvKEKaXOUH4Cdq6AbQXebW8h3PB/0Pta28mUSgpa6FXicSOwZ220sL8BHy2HyClw0qDzYGjdDf5xB2RdBOd1t51WqYSnhV4lhnAZvPcX2LoEdrwFJw97z2f1hSG3QbfL4YI8aNQMDn0If7gM/jYFvrkI0htbja5UotNCrxLDmw/Am/8LLbtC7+uge75X3JtVMUdT6wvgC4/DrMmw6Ecw8cEGj6tUMtFCr+wr3Q/Lfg99vgBf/hOI1LxO7gTI+y4s/a23p9/3+ninVCpp6agbZd/bD0H4BFzxI39F/rSrfgpdhsKc78GB4vjlUyrJaaFXdh0ugZVPwcU3QduetVs3lO6NvgllwN9u9UbnKKXOooVe2fXGA2BcuPwHdVu/ZWf40kzYux5eqeN7KBVwWuiVPZ9shfeegUHf8A6w1lXOGBh5D6x5Gt6fFbt8SgWEFnplT8H/eN0ul91b//e64kdwwQh4+fuwb2P930+pANFCr+zYWwjrXoCh34bmWfV/v1AaXP8UpGd6/fVlx+r/nkoFhBZ6ZceSX0Cj5jDirti9Z4sOcP2TsH8TzJuuc+IoFaWFXjW8ktWwaZ43Dj7zvNi+94VXwOU/hPef9/r/lVJa6JUFS34OmefDsO/E5/0v/4F3Vu38e2FfUXw+Q6kkooVeNaztb3mTlY28x+u6iQcn5HXhOOnw9sPx+QylkogWetVwjIEl/wXNO8Dgf4nvZzVrB/2/AoUvwfGD8f0spRKcFnrVcLYs8uaVv/wHkN4k/p936RRveuMPZsf/s5RKYFroVcNwXW9vvnU2DPh6w3xmh/7QcSCs/pOOwFEpTQu9ahgb/gEfr4P8f/fmqGkol06B/UWw892G+0ylEoxOU6ziTtwIvP7f0DYX+t1Qr/fKnjGvVstn0ox3GzVmwRP3cW/5NN/r7bh/Ym2jKZWwdI9exV3W3tfhky1w5Y+9ETEN6DiNmRPJY6KznBbo2bIqNWmhV/EVPkX2jr9CxwGQe42VCM9FrqSJlDEp9I6Vz1fKNi30Kr7WPkfjU/vgyp/U7qIiMbTedGedm81NodcAPSirUo+vQi8i40Rkk4gUi8iMKl7PFZFlInJKRO6tzboq4Nb/nWOZneHCK63GeD5yFb2dnVwiW63mUMqGGgu9iISAR4HxQB9gsoj0qbTYQeB7wK/qsK4KquMH4cOlHGgzzNre/GlzIsM5ZhoxObTEag6lbPCzRz8EKDbGbDPGlAGzgEkVFzDG7DPGrATKa7uuCrAti8FEONBmiO0klJLJnEge14aW0ZzjtuMo1aD8FPpOwM4Kj0uiz/lRn3VVstv4MjRrz9HmObaTAPB85Eoy5ZQelFUpx884+qr+5vZ7RMv3uiIyFZgKkJWVRUFBgc+POFNpaWmd1000ydwWJ1LGiE2L2Jt1OaXHjsesHdP7heu+sunK/p1duavRa2R1yT9nd1J1eZN5m1QWlLYEpR0Qv7b4KfQlQJcKjzsDu32+v+91jTEzgZkAgwYNMvn5+T4/4kwFBQXUdd1Ek9Rt2bwI3JN0vPI2Nu9Kj1k7ptTyhKnK9oau4r70P7KocCfrTPdql9vxtfwqn0/qbVJJUNoSlHZA/Nrip+tmJZAjIt1EJAO4EZjj8/3rs65KZpvmQUYz6DbKdpIz/DMyghMmg8mh12xHUarB1FjojTFh4E5gIVAEzDbGFIrINBGZBiAi7UWkBLgH+LGIlIhIi+rWjVdjVIJwXdj0CvS4CtIa2U5zhqNkMjcynEmhpTTlhO04SjUIX3PdGGPmA/MrPfd4hfsf43XL+FpXBdzuNVC6F3ol5nwxz0eu5Ctpb3BdaCnPR66yHUepuNMzY1XsbXwZJAQ9x9pOUqX3TA+K3C46pl6lDC30KvY2zofsEdCkte0k1RCej1xJf2c7fWWb7TBKxZ0WehVbn2yFA5sSttvmtH9ERnLSpDM59LrtKErFnRZ6FVsbo8MfcyfYzVGDIzTlZXc414WWkslJ23GUiist9Cq2Ns2HrH7QqqvtJDV6PnwFzeUE14SW2Y6iVFxpoVexc+yAd/HvBN+bP2216clmtxM36UFZFXB6KcEAq+1l9+rqs8vubV4AxoVeyVHoTx+U/Wn6X+gjO9hgsm0HUioudI9exc7GedCiM3S42HYS316MXMYpk86XQ2/YjqJU3GihV7FRdhy2vu5121iee742DtOMN9z+jAutRHBtx1EqLrTQq9jY9jqETyRRt83nFkQG00EO0l/H1KuA0kKvYmPjfGjUErJH2k5Sa6+5Awkbh6tDq2xHUSoutNCr+nMj3oHYnDEQSredptYO04xlbh+udlaiFw9XQaSFXtXfznfh+IGkGVZZlYXuYC509tBDdtmOolTMaaFX9bfxZXDSoccY20nqbFFkEADjnJWWkygVe1roVT0Z72zYbqOgcQvbYepsH61Z4/bg6pAWehU8WuhVvfSQXXBwW1J325y2IDKYfs4OOst+21GUiikt9KpexjqrvTtJOKyysoXuYADGOjr6RgWLFnpVL2NCq6HjAGjR0XaUevvQtKfI7aLdNypwtNCrOmvLIQY4xQk/93xtLHIHM1g2Qek+21GUihkt9KrORofWeHcC0D9/2oLIYByJHmBWKiC00Ks6G+us4kO3HbTrYztKzBSZrnzktoWil21HUSpmtNCrOmnKCfKcQha7lybVJGY1Exa4Q2BbAZw8bDuMUjGhhV7VySjnAxpJmMXRE42CZGFkELjlsGWx7ShKxYQWelUno0NrOGSascr0tB0l5taYHGiWBUVzbEdRKia00Ks6MOQ5hbzj9iVCyHaYmDM43nkBW16F8hO24yhVb1roVa11k4/pIAdZ6l5kO0r89L4Wyo95F1NRKslpoVe1NtzZAMAyNzijbc6SfZk3v/5GHX2jkp8WelVreU4he8x5bDftbUeJn7QM6DXOG08fCdtOo1S9+Cr0IjJORDaJSLGIzKjidRGR30Rf/0BEBlZ47fsiUigi60XkeRFpHMsGqIYluAxzNrDU7QMEaVhlFXKvgROH4MN3bCdRql5qLPQiEgIeBcYDfYDJIlL5b/bxQE70NhV4LLpuJ+B7wCBjTF8gBNwYs/SqwfWUEtrIEZYFuX/+tB5XQVoTKJprO4lS9eJnj34IUGyM2WaMKQNmAZMqLTMJeNp4lgOtRKRD9LU0oImIpAGZwO4YZVcWfNY/Hwlw//xpGU29Yr9xHriu7TRK1Vmaj2U6ATsrPC4BhvpYppMxZpWI/Ar4CDgBLDLGLKrqQ0RkKt5fA2RlZVFQUOCrAZWVlpbWed1EU9+2TO8X+77libvXc7isHTfmtAa8968pYyy3STzaVJXTebPoQe+jL7P65Sc42qKX/n4loKC0A+LXFj+FvqqO2MpXUK5yGRFpjbe33w34FPibiNxsjHnmrIWNmQnMBBg0aJDJz8/3Ee1sBQUF1HXdRFPftkyZMS92YQAHl1sbbWReZCgPrvv8V2fH1/LPuV4st0ms21Sdz9p04mLY/DsubbIb8r+tv18JKCjtgPi1xU/XTQnQpcLjzpzd/VLdMqOB7caY/caYcuBFIK/ucZVNF8kOWsjx1OifP61Ja2+oZdFcMJX3b5RKDn4K/UogR0S6iUgG3sHUyueGzwFuiY6+GQYcNsbsweuyGSYimSIiwFVAUQzzqwY03CkEYJnb23KSBtb7Gji4FfZvtJ1EqTqpsdAbY8LAncBCvCI92xhTKCLTRGRadLH5wDagGHgCuD267grgBWANsC76eTNj3QjVMPKcDWxxO7Gf1rajNKzTF1bR0TcqSfnpo8cYMx+vmFd87vEK9w1wRzXr/hT4aT0yqgSQRpjBzkZeiIyyHaXhtegAnYd4hT53iO00StWanhmrfOkv22gqp1Krf76i3tfAxx/Q+MRe20mUqjUt9MqXvGj//PJU658/LfcaANocWG45iFK1p4Ve+ZLnFFLoXsCnNLcdxY7zL4R2F2mhV0lJC72qUSPKuNTZEuxpif3InUjLwxvh2AHbSZSqFS30qkYDnS00kvJgT0vsR+5EBBc2L7CdRKla0UKvajTcKSRsHN51c21HsavDxZxs1AY2zq95WaUSiBZ6VaM8ZwPrTHdKybQdxS4RPjl/CGxdAmXHbadRyjct9OqcMjnJxbI1Ov+8OtBmKIRPwDa9xKBKHlro1TkNdjaRLpHUHT9fyaet+kYvMajdNyp5aKFX5zTcKaTMhFjl9rQdJSEYJw1yxsDmV8CN2I6jlC9a6NU55TmFvGdyOEkj21ESR+5EOP4J7FxhO4lSvmihV9VqQSl9ZYcOq6ysx2hw0r0rTymVBLTQq2oNdTbiiGFpRPvnz9C4BXQb5RV6naNeJQEt9KpaeU4hJ0wGa00P21EST+5EOLRd56hXSUELvarWcGcDK91elJFuO0ri6TXB+6ndNyoJaKFXVTqfw+Q6O1mu/fNVa9EBOg6ETTrMUiU+LfSqSsMc74qPKT+R2bnkToRdq+HIHttJlDonLfSqSnlOIUdNE9aZbrajJK7c6CUGda9eJTgt9KpKw51CVri5RAjZjpK42uZC625a6FXC83XNWJVa2vMJ3Z2Pebb8Kl/LZ8849wHJ6f3CTKlhmaQk4u3Vr/gDnDziDbtUKgHpHr06y3BnA4DOb+NH7kRwy6H4VdtJlKqWFnp1ljynkEOmGUWmq+0oia/LUMg8X7tvVELTQq8qMQwPbWC52xujvx41c0LQczxsXgSRcttplKqSfpPVGbrIPjrLAR1WWRu5E+DUYdjxtu0kSlVJC706Q160f14LfS10vwLSmmj3jUpYWujVGfKcQvaZVmw1HW1HSR4ZmXDhld7FSHSSM5WAtNCrCgx5TmF0WmKxHSa55E6AIyWw533bSZQ6i69CLyLjRGSTiBSLyIwqXhcR+U309Q9EZGCF11qJyAsislFEikRkeCwboGKnp5TQVg7zttvXdpTk03MciKPdNyoh1VjoRSQEPAqMB/oAk0Wk8kxX44Gc6G0q8FiF1x4BFhhjcoGLgaIY5FZxMMJZD6Dzz9dF0zbQZZjOZqkSkp89+iFAsTFmmzGmDJgFTKq0zCTgaeNZDrQSkQ4i0gIYBTwFYIwpM8Z8Grv4KpZGOOvZ7maxi7a2oySn3Amwdz0c2mE7iVJn8DMFQidgZ4XHJcBQH8t0AsLAfuCPInIxsBq4yxhzrPKHiMhUvL8GyMrKoqCgwGcTzlRaWlrndRNNfdsyvV/Y97KOCTNqWxGbmg1nepb/9fzIalK7LImgun/3c22TJsfbMBTYMu+37Op8bdyyxUpQvitBaQfEry1+Cn1VR+UqDy2obpk0YCDwXWPMChF5BJgB/OSshY2ZCcwEGDRokMnPz/cR7WwFBQXUdd1EU9+21GZ+mYGyjTsaneT3B/rzyr7YToE0vV+YB9cl17RKO76WX+XzNW6T7Y+QE95ETv6DcckVS0H5rgSlHRC/tvjpuikBulR43BnY7XOZEqDEGLMi+vwLeIVfJZgRznpcI3oh8PrKnQAfLoXjB20nUeozfgr9SiBHRLqJSAZwIzCn0jJzgFuio2+GAYeNMXuMMR8DO0WkV3S5q4ANsQqvYmdkaD2F5gI+pbntKMktdyKYCGxZZDuJUp+psdAbY8LAncBCvBEzs40xhSIyTUSmRRebD2wDioEngNsrvMV3gWdF5APgEuC/YxdfxUITTjJAtvCODqusvw4DoHkH2Piy7SRKfcZXx6kxZj5eMa/43OMV7hvgjmrWXQsMqntEFW9DnE1kSESnPYgFx/H26t97Fk6VQqNmthMppWfGKshz1nPKpLHS7VXzwqpmF30Jwidg8wLbSZQCtNArYIRTyHsmhxM0th0lGLoO97pv1r9oO4lSgBb6lNeaI/R1dvB2RPvnY8ZxvL364sVw4lPbaZTSQp/qhuu0xPHR93qIlOmUCCohaKFPcSOd9Rw1TXjfXGg7SrB0GgitLoD1f7edRCkt9KkuzylkudubCCHbUYJFxNur31YAxw7YTqNSnBb6FNZZ9pPt7NXx8/HS93rv5KkN/7SdRKU4LfQpLC86LbEW+jjJugja9NLRN8o6LfQpbKSznn2mFVtMJ9tRgul0982H78CRytNDKdVwtNCnLMNwp5B33IvQywbGUd/rAQOFL9lOolKYFvoU1Ut20laO6LDKeGvTAzpcrKNvlFVa6FPUyNP983qiVPz1vR52rYaD220nUSlKC32KynMK2ea2ZzdtbEcJvou+6P0s1IOyyo7kuuyPiok0wgx1ingpMtJ2lISVXc3Vuab3C9fqyl2nvZDRk6aL/8T4eblnPL/j/ol1yqdUbegefQq6WLbSTE7qsMoGNDcynN7OR/SQEttRVArSQp+CRuplAxvc/MhQIka4NrTcdhSVgrTQp6C8UCHrTTaH0YtiNJT9tGK524drnGWAsR1HpRgt9CkmM3rZwKXabdPg5rrDudDZw0Xyoe0oKsVooU8xQ5yNZEiEt7XQN7gFkcGUmxDXhpbZjqJSjBb6FJPnFHLKpLHK7Wk7Ssr5lOa85fbjmpB236iGpYU+xYx01rPa7clJGtmOkpLmRobTWQ4wULbYjqJSiBb6FHIeR+jjfKjDKi1a7F7KKZOu3TeqQWmhTyF5TiGg0xLbVEomS9xLmBhagYNrO45KEVroU0ies54jpgnrTDfbUVLa3Mhw2smnDHWKbEdRKUILfQoZ6axnudtHLxto2RJ3AKWmMdc6S21HUSlCC32K6Cz76Ors126bBHCSRix2L2V8aCWEy2zHUSlAC32KGPFZ/7zOP58I5kaG01pKvYuHKxVnvgq9iIwTkU0iUiwiM6p4XUTkN9HXPxCRgZVeD4nIeyLycqyCq9oZ6axjr2lFsV42MCG85fbnsMnUC5KoBlFjoReREPAoMB7oA0wWkcqzYY0HcqK3qcBjlV6/C9AjT5Y4uOQ5hdFuG71sYCIoJ41XIkNg4zwoO2Y7jgo4P3v0Q4BiY8w2Y0wZMAuYVGmZScDTxrMcaCUiHQBEpDMwEXgyhrlVLQySTZwvR3ktMrDmhVWDmR3Jh7Kj8MFs21FUwPm58EgnYGeFxyXAUB/LdAL2AA8DPwCan+tDRGQq3l8DZGVlUVBQ4CPa2UpLS+u8bqKpb1um9wsDcNn+FYQPp9Ovd19ynXCM0vmX1eTzLMkupm0x3Tj6cXfk9V+z6mg2SMP+tRWU70pQ2gHxa4ufQl/Vb1/liTqqXEZErgH2GWNWi0j+uT7EGDMTmAkwaNAgk59/zsWrVVBQQF3XTTT1bcuUGfMQXK5vtIrX3f7cX2hnWuLp/cI8uC4YFzOLdVu++5V7YM6d5HdLh+yGveJXUL4rQWkHxK8tfrpuSoAuFR53Bnb7XGYEcJ2I7MDr8rlSRJ6pc1pVa/1lGx3lIK9EBtuOoqrS7wZo0hrenWk7iQowP7smK4EcEekG7AJuBG6qtMwc4E4RmYXXrXPYGLMH+Lfojege/b3GmJtjE135MT60knIT4lVX++cTUfZPljAjbQTfKpzLyBlP8zHnx+2z9Pq0qavGPXpjTBi4E1iIN3JmtjGmUESmici06GLzgW1AMfAEcHuc8qpaMVztvMsytw9H9GpSCeuZyGgcDDelvWY7igooX52Nxpj5eMW84nOPV7hvgDtqeI8CoKDWCVWd5cpOujl7mVl+je0o6hxKTDtecwcwObSE34W/SBnptiOpgNEzYwNsfOhdXCMsjgyyHUXV4OnIWNrKEcY7K2xHUQGkhT7ArnZWstL04gAtbUdRNXjb7ctWtwO3pi2yHUUFkBb6oDpQTK6zkwU62iYpGBz+EhnDQKeYfrLNdhwVMFrog6poDgALIkMsB1F+/T0yimOmEbeEdK9exZYW+qAqmsNa90L2xHG4noqto2TyYuQyrgstozVHbMdRAaKFPog+/Qh2v6fdNkno6chYGkk5Xw0V2I6iAkQLfRAVebNBv+Jqt02y2WI6szTSh5vTXtVryqqY0UIfREVzIKsvH5r2tpOoOvhzZCyd5QBXOWtsR1EBoYU+aI7uhY+WQ+9rbSdRdfSqeym7zXl6UFbFjBb6oNn4MmCg93W2k6g6ihDimfBoLgut50LZZTuOCgAt9EFTNAfOuxDa9badRNXDXyNXcMqk8fXQYttRVABooQ+S4wdh+1vQ57oGv4iFiq1PaMnL7jBuCL1JM47bjqOSnBb6INn0CpiIdtsExNPhsTSTk3wx9LbtKCrJaaEPkqI50LILdBxgO4mKgfdND9a63bk1tIizL+qmlH9a6IPi1FHYusQbbaPdNoHxl/BYeji7yXMKbUdRSUwLfVBsXgiRMu22CZiX3WEcNM2ie/VK1Y0W+qAomgtN20EXPRs2SE6RwazIlYx2VpMte2zHUUlKC30QlJ+ALYuh9zXghGynUTH2x/A4TpLBD9Nm2Y6ikpQW+iAofg3Kj+nZsAG1n1Y8Fr6O8aGVDJEi23FUEtJCHwRFc6FxK8i+zHYSFSdPRiaw25zHj9OfQXSyM1VLWuiTXbjMGz+fOxFCelHpoDpJI35ZfiP9ne180dFx9ap2tNAnu+1vwqnD2m2TAua4eax1u/Ov6bNpwknbcVQS0UKf7IrmQEYz6H6F7SQqzgwO95XfTAc5yG2h+bbjqCSihT6ZRcKwcR70vBrSG9tOoxrAKpPLvMgQpqXNpR2HbMdRSUILfTIr+iccPwB9r7edRDWgX4YnEyLCvWmzbUdRSUILfbIyBt76NZyfAz3H206jGtBHJos/Ra7mhtCbXCTbbcdRSUALfbIqfg32roORd4OjmzHVPBr+Aodoxo/TnkUnPFM1SfOzkIiMAx4BQsCTxpj7K70u0dcnAMeBKcaYNSLSBXgaaA+4wExjzCMxzJ90smfM873s9H5hplSz/KyM/6KrnMflf21O+V/9v6cKhiM05dfhG7gv/Y+MiaxmsTvIdiSVwGrcFRSREPAoMB7oA0wWkT6VFhsP5ERvU4HHos+HgenGmN7AMOCOKtZVtTRQNjPMKeLJ8ETK/f1frQLo+ciVbHE78W9pz5FO2HYclcD8/M0/BCg2xmwzxpQBs4BJlZaZBDxtPMuBViLSwRizxxizBsAYcxQoAjrFMH9KmpY2l0OmGbMiOqQylUUI8YvwTXR3PuZmveSgOgc/hb4TsLPC4xLOLtY1LiMi2cAAYEWtU6rP9JASxoZW8+fIWI6jQypTXYF7CW9G+nFX2ou0pNR2HJWg/PzdX9VVLCof/TnnMiLSDPg7cLcx5kiVHyIyFa/bh6ysLAoKCnxEO1tpaWmd120I0/v5/xM7q8nZy4/+eA7lpRk07XkV00PJ8ed6Ve1IVonYlpJTNzLyox/z5/Yv8Fbbm6tdrvL3ItG/K34FpR0Qv7b4KfQlQJcKjzsDu/0uIyLpeEX+WWPMi9V9iDFmJjATYNCgQSY/P99HtLMVFBRQ13UbQnUHV6syvV+YB9d9vok6coBpjZbxl8gYfrGhdTzixUXldiSzxGxLN0jL58uHXuP7e69mu+lQ5VI7vpZ/xuNE/674FZR2QPza4qfrZiWQIyLdRCQDuBGYU2mZOcAt4hkGHDbG7ImOxnkKKDLGPBTT5CnotjTvP4knwxMsJ1GJ5qHwVzhFOv+W9pztKCoB1VjojTFh4E5gId7B1NnGmEIRmSYi06KLzQe2AcXAE8Dt0edHAF8HrhSRtdGbVqk6aM0Rbgy9zj8iI9hNG9txVII5QEt+H57E2NBqRjurbcdRCcbX36DGmPl4xbzic49XuG+AO6pY722q7r9XtTQlbRFNpIzHIzpLparaU5HxTAit4KH03/PFsp+z1egAN+XRUyqTQCYnuTW0kIWRQfrlVdU6RQbfLvs+5aQxM/0hmnPcdiSVILTQJ4HJoSW0kmM8Hta9eXVuu2jL7WV301X28ev0R/VqVArQQp/w0gnzrbT5LIv04T2TYzuOSgIrTG9+Hv46o0PvcU/aC7bjqASghT7BfSH0Nh3kII9p37yqhb9ExjArnM930/7BeEfPUUx1iTYgWFUgxmVaaC7r3WzedPvbjqOSivAf4W+Q4+ziwfTH2V5W9dh6lRp0jz6BdT+2mgudPTwWvg4dvKRqq4x0ppXdzREymZn+IBw/aDuSskQLfcIyXHpwHjvcLF5xh9gOo5LUflozrez7ZMkh+Nut3uUnVcrRQp+ghjsbyDq1jT9ErsHVzaTqYa3pwY/C/wLb34TFP7EdR1mgFSRBfSc0h2OhlrwYucx2FBUAL0Quh6HTYPnvYe3ztuOoBqYHYxPQKOd9RoXW8U6rr3LqWIbtOCooxt4Hewth7l3QtqftNKoB6R59gsmWPfw2/bcUuV34oNVo23FUkITS4ct/huZZMOtmMk59YjuRaiC6R59AmnOcJ9IfIoLDbeXT+arTyHYkFTRNz4cbn4Onrmbgmn+Fvt2g00DbqRJeba71XB9/Gtc0Lu+re/QJwsHl4fRHyZaPub38bkpMO9uRVFC17wfffAVw4I/j4f2/2k6k4kwLfYK4N202V4Xe4z/Dt7Dc1eunqzjrcDGrL30QOg+Gl6bCwh/p0MsA00KfAK5z3uH2tDk8G76KZyLaL68aRnlGS/j6SzDk27Dsd/Ds9XpSVUBpobesn2zjgfSZrHBz+Vn4VvQMWNWgQukw4QGY9Ch8uBRm5nsjc1SgaKG3qC2HmJnxEAdoyXfK7qZcj40rWwbcDFPmQ/gUPDkGNvzTdiIVQ1roLcmgnD9k/JqWHOO2sukcpIXtSCrVdRkM334DsvrA7Fvgtf8CV+ezDwIt9FYYfpH2FAOdYu4p/w5F5gLbgZTyNG8PU+bBgK/DW7+CWZPhxKe2U6l60kJvwTdDC/hy2ps8Ev4SC3TCMpVo0hrBdb+FCb+C4lfhN5fAm/8LJ4/YTqbqSAt9A7vM+YAfpT3DgshgHg5/yXYcpaomAkNug2+9Bl2GwpL74OF+8MYDcPKw7XSqlrTQN6QPl/G79N+w2XTmnvLvYPSfXyW6jpfATX+F216HrsPh9V94Bb/gfu3SSSJaaRrCpx/B374BfxxHKU24rXw6x2lsO5VS/nUaCDfNgqlvwAUjoeB/4OH+8Pr/wIlDttOpGuh4vng6VQrvPAxLfwsIXD6D0Qt7ckKLvEpWHS+Byc/Bng/gjV/CG/d7Ux8Pnebdmp5vO6Gqghb6eHBdWDcbXv0ZHN0DfW+AMf8JLTtzYmHDTI6kVFx16A83Pgsfr/MK/psPeKN0Og+GnDGQMxba9/f6+pV1Wuhjbee7sGAG7FoNHQd608J2HWo7lVLx0b4ffPUZ2FcEhS/BlkXegdsl90Gz9pAz2iv63fOhcUvbaVOWFvpYOVzi7cGv+5v3C/6Fx6H/V8HRwyAqBbTr7d2u+Hco3ecNy9yyCDbMhfeeASfNO5ibMwa6jYI2vSAj03bqlKGFvq7KjsGe96Fklbf3vnkhYGDUv8KIu6FRM9sJlbKjWTu45CbvFglDybte0d+yGBb/x+fLtezqXemqTa8KP3tB5nn2sgeUr0IvIuOAR4AQ8KQx5v5Kr0v09QnAcWCKMWaNn3WTQiQM+4u8gr5rNexaA/s2gImeHt6qK/S7Hi7/oXdfKeUJpcEFed5t9M/g8C4oWQkHNsP+TXBgE+x4B8InPl8ns41X8M+/EJq2g6ZtoGlbyDzf+3n6fkj3U/2q8V9KRELAo8AYoARYKSJzjDEbKiw2HsiJ3oYCjwFDfa4bP8Z4xdiNgImceb/8BJw66p3td+pw9OeRCs8d8X4e2g67137+i9ikNXS6FHInej87DoRmbRukOUolvZadvFtFrguHP4L9m73Cv3+T9x/BpgVw/BPv+1qVJq0hsw0DytNgZ0evKyi9afRnJmQ0hfQmZz4XyvBm7HTSKtxP9/7TcNKjz6WBhMAJgTggIVpzBBcnehMiOBgEFwcDGCR68+4n2iy0fv5LHAIUG2O2AYjILGASULFYTwKeNsYYYLmItBKRDkC2j3Vj55fdGHXyCLyJV9AxdXwjgUYtoHELaN4BLp3iFfVOA+G87jqSQKlYchxone3deo498zXXhZOfwrEDcGw/HI/+PPZJ9Od+3N1bvWWO7IbyY1B2HMqjtxh5r44joiPG+w9glenFjWU/iVme2vJT6DsBOys8LsHba69pmU4+1wVARKYCU6MPS0Vkk49sVWkDHKjjuhWcPs27EHi1/m9XB9+LWVvsCko7ILnbIr8866mkbUslSdCOFcA1NS51xS/r1ZZqZ0f0U+ir2n2tvKtc3TJ+1vWeNGYmMNNHnnMSkVXGmEH1fZ9EEJS2BKUdoG1JREFpB8SvLX4KfQnQpcLjzsBun8tk+FhXKaVUHPkZ5L0SyBGRbiKSAdwIzKm0zBzgFvEMAw4bY/b4XFcppVQc1bhHb4wJi8idwEK8IZL/Z4wpFJFp0dcfB+bjDa0sxhte+Y1zrRuXlnyu3t0/CSQobQlKO0DbkoiC0g6IU1vEGyijlFIqqPT8fKWUCjgt9EopFXCBKfQiMk5ENolIsYjMsJ2nPkRkh4isE5G1IrLKdp7aEJH/E5F9IrK+wnPnichiEdkS/dnaZka/qmnLz0RkV3TbrBWRCTYz+iEiXUTkdREpEpFCEbkr+nzSbZdztCUZt0tjEXlXRN6PtuU/o8/HfLsEoo8+OtXCZipMtQBMbrCpFmJMRHYAg4wxCX4SyNlEZBRQinemdN/ocw8AB40x90f/E25tjPmhzZx+VNOWnwGlxphf2cxWG9Gz1DsYY9aISHNgNfAFYApJtl3O0ZavkHzbRYCmxphSEUkH3gbuAr5EjLdLUPboP5umwRhTBpyeakE1MGPMm8DBSk9PAv4cvf9nvC9mwqumLUnHGLPn9CSDxpijQBHeWetJt13O0ZakYzyl0Yfp0ZshDtslKIW+uikYkpUBFonI6ujUEMkuK3peBdGf7Sznqa87ReSDaNdOwnd3VCQi2cAAvHPyk3q7VGoLJOF2EZGQiKwF9gGLjTFx2S5BKfS+p1pIEiOMMQPxZgW9I9qFoBLDY8CFwCXAHuBBq2lqQUSaAX8H7jbGHLGdpz6qaEtSbhdjTMQYcwnerAFDRKRvPD4nKIXezzQNScMYszv6cx/wEl7XVDLbG+1bPd3Hus9ynjozxuyNfjld4AmSZNtE+4D/DjxrjHkx+nRSbpeq2pKs2+U0Y8ynQAEwjjhsl6AU+sBMtSAiTaMHmRCRpsBYYP2510p4c4Bbo/dvBf5pMUu9nP4CRn2RJNg20YN+TwFFxpiHKryUdNulurYk6XZpKyKtovebAKOBjcRhuwRi1A1AdDjVw3w+1cIv7CaqGxHpjrcXD94UFc8lU1tE5HkgH2/q2L3AT4F/ALOBrsBHwJeNMQl/kLOatuTjdQ8YYAfw7dP9qYlKREYCbwHrgOhl0fh3vL7tpNou52jLZJJvu/THO9gawtvpnm2M+bmInE+Mt0tgCr1SSqmqBaXrRimlVDW00CulVMBpoVdKqYDTQq+UUgGnhV4ppQJOC71SSgWcFnqllAq4/wfIiMA2NShecwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "k = np.arange(0,30)\n",
    "mu = np.exp(result.params.Intercept)\n",
    "pmf_poisson = stats.poisson.pmf(k,mu)\n",
    "\n",
    "plt.figure()\n",
    "df['NamedStorms'].hist(density=True)\n",
    "plt.plot(k,pmf_poisson)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Poisson regression model, including a linear temporal trend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Generalized Linear Model Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>      <td>NamedStorms</td>   <th>  No. Observations:  </th>  <td>    70</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                  <td>GLM</td>       <th>  Df Residuals:      </th>  <td>    65</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model Family:</th>         <td>Poisson</td>     <th>  Df Model:          </th>  <td>     4</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Link Function:</th>          <td>log</td>       <th>  Scale:             </th> <td>  1.0000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -186.77</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Sun, 14 Mar 2021</td> <th>  Deviance:          </th> <td>  73.614</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>21:48:33</td>     <th>  Pearson chi2:      </th>  <td>  73.3</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>          <td>4</td>        <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>     <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   -8.9786</td> <td>    3.621</td> <td>   -2.480</td> <td> 0.013</td> <td>  -16.076</td> <td>   -1.881</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Year</th>      <td>    0.0057</td> <td>    0.002</td> <td>    3.154</td> <td> 0.002</td> <td>    0.002</td> <td>    0.009</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>NAO</th>       <td>   -0.0267</td> <td>    0.048</td> <td>   -0.558</td> <td> 0.577</td> <td>   -0.120</td> <td>    0.067</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>SOI</th>       <td>    0.0620</td> <td>    0.037</td> <td>    1.654</td> <td> 0.098</td> <td>   -0.011</td> <td>    0.135</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>TNA</th>       <td>    0.3563</td> <td>    0.096</td> <td>    3.729</td> <td> 0.000</td> <td>    0.169</td> <td>    0.544</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                 Generalized Linear Model Regression Results                  \n",
       "==============================================================================\n",
       "Dep. Variable:            NamedStorms   No. Observations:                   70\n",
       "Model:                            GLM   Df Residuals:                       65\n",
       "Model Family:                 Poisson   Df Model:                            4\n",
       "Link Function:                    log   Scale:                          1.0000\n",
       "Method:                          IRLS   Log-Likelihood:                -186.77\n",
       "Date:                Sun, 14 Mar 2021   Deviance:                       73.614\n",
       "Time:                        21:48:33   Pearson chi2:                     73.3\n",
       "No. Iterations:                     4                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -8.9786      3.621     -2.480      0.013     -16.076      -1.881\n",
       "Year           0.0057      0.002      3.154      0.002       0.002       0.009\n",
       "NAO           -0.0267      0.048     -0.558      0.577      -0.120       0.067\n",
       "SOI            0.0620      0.037      1.654      0.098      -0.011       0.135\n",
       "TNA            0.3563      0.096      3.729      0.000       0.169       0.544\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = smf.glm(formula='NamedStorms ~ Year + NAO + SOI + TNA',\n",
    "                 data=df,\n",
    "                 family=sm.families.Poisson()).fit()\n",
    "result.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta = result.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept   -8.978558\n",
       "Year         0.005747\n",
       "NAO         -0.026677\n",
       "SOI          0.061967\n",
       "TNA          0.356328\n",
       "dtype: float64"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercises\n",
    "* Create a Python variable `y` that contains the observed count data.\n",
    "* Create a Python variable `yhat` that contains the modeled count data (use the parameters in `beta`).\n",
    "* Make plots comparing the observed and modeled counts."
   ]
  }
 ],
 "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}