maximum_likelihood_GLM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Maximum likelihood & the GLM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "import numpy as np\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th> <td>   0.100</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>  -0.002</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>  0.9843</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 27 Jun 2018</td> <th>  Prob (F-statistic):</th>  <td> 0.463</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>16:27:10</td>     <th>  Log-Likelihood:    </th> <td> -120.09</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   100</td>      <th>  AIC:               </th> <td>   262.2</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    89</td>      <th>  BIC:               </th> <td>   290.8</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>    10</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>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>    0.0736</td> <td>    0.092</td> <td>    0.801</td> <td> 0.425</td> <td>   -0.109</td> <td>    0.256</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>    <td>    0.0218</td> <td>    0.084</td> <td>    0.260</td> <td> 0.795</td> <td>   -0.145</td> <td>    0.189</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x2</th>    <td>   -0.1252</td> <td>    0.093</td> <td>   -1.350</td> <td> 0.180</td> <td>   -0.310</td> <td>    0.059</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x3</th>    <td>   -0.0367</td> <td>    0.081</td> <td>   -0.455</td> <td> 0.650</td> <td>   -0.197</td> <td>    0.124</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x4</th>    <td>    0.0511</td> <td>    0.088</td> <td>    0.580</td> <td> 0.564</td> <td>   -0.124</td> <td>    0.227</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x5</th>    <td>    0.0819</td> <td>    0.087</td> <td>    0.938</td> <td> 0.351</td> <td>   -0.092</td> <td>    0.255</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x6</th>    <td>   -0.0537</td> <td>    0.084</td> <td>   -0.636</td> <td> 0.526</td> <td>   -0.221</td> <td>    0.114</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x7</th>    <td>   -0.1436</td> <td>    0.087</td> <td>   -1.660</td> <td> 0.101</td> <td>   -0.316</td> <td>    0.028</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x8</th>    <td>   -0.1349</td> <td>    0.089</td> <td>   -1.517</td> <td> 0.133</td> <td>   -0.312</td> <td>    0.042</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x9</th>    <td>    0.0075</td> <td>    0.079</td> <td>    0.095</td> <td> 0.925</td> <td>   -0.149</td> <td>    0.164</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x10</th>   <td>    0.1311</td> <td>    0.100</td> <td>    1.306</td> <td> 0.195</td> <td>   -0.068</td> <td>    0.331</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 0.453</td> <th>  Durbin-Watson:     </th> <td>   1.995</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.798</td> <th>  Jarque-Bera (JB):  </th> <td>   0.559</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.149</td> <th>  Prob(JB):          </th> <td>   0.756</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.787</td> <th>  Cond. No.          </th> <td>    1.92</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   R-squared:                       0.100\n",
       "Model:                            OLS   Adj. R-squared:                 -0.002\n",
       "Method:                 Least Squares   F-statistic:                    0.9843\n",
       "Date:                Wed, 27 Jun 2018   Prob (F-statistic):              0.463\n",
       "Time:                        16:27:10   Log-Likelihood:                -120.09\n",
       "No. Observations:                 100   AIC:                             262.2\n",
       "Df Residuals:                      89   BIC:                             290.8\n",
       "Df Model:                          10                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          0.0736      0.092      0.801      0.425      -0.109       0.256\n",
       "x1             0.0218      0.084      0.260      0.795      -0.145       0.189\n",
       "x2            -0.1252      0.093     -1.350      0.180      -0.310       0.059\n",
       "x3            -0.0367      0.081     -0.455      0.650      -0.197       0.124\n",
       "x4             0.0511      0.088      0.580      0.564      -0.124       0.227\n",
       "x5             0.0819      0.087      0.938      0.351      -0.092       0.255\n",
       "x6            -0.0537      0.084     -0.636      0.526      -0.221       0.114\n",
       "x7            -0.1436      0.087     -1.660      0.101      -0.316       0.028\n",
       "x8            -0.1349      0.089     -1.517      0.133      -0.312       0.042\n",
       "x9             0.0075      0.079      0.095      0.925      -0.149       0.164\n",
       "x10            0.1311      0.100      1.306      0.195      -0.068       0.331\n",
       "==============================================================================\n",
       "Omnibus:                        0.453   Durbin-Watson:                   1.995\n",
       "Prob(Omnibus):                  0.798   Jarque-Bera (JB):                0.559\n",
       "Skew:                           0.149   Prob(JB):                        0.756\n",
       "Kurtosis:                       2.787   Cond. No.                         1.92\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 100\n",
    "K = 10\n",
    "X = np.random.normal(0, 1, (N, K))\n",
    "X = np.c_[np.ones(N), X]\n",
    "y = np.random.normal(0, 1, N)\n",
    "\n",
    "model = sm.OLS(y, X)\n",
    "fitted_mod = model.fit()\n",
    "mod.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-139.37828006536603\n",
      "-139.37828006536603\n"
     ]
    }
   ],
   "source": [
    "betas = np.array(mod.params)\n",
    "sigma_sq = np.sqrt(mod.resid.dot(mod.resid) / (N))# - K))\n",
    "llh = np.log(np.prod(scipy.stats.norm(X.dot(betas), sigma_sq).pdf(y)))\n",
    "print(llh)\n",
    "llh = scipy.stats.norm(X.dot(betas), sigma_sq).logpdf(y).sum()\n",
    "print(llh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "30"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t = np.array([3, 4, 1, 2])\n",
    "t.dot(t.T)"
   ]
  }
 ],
 "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
}