Chapter-13/13.3 Introduction to statsmodels（statsmodels简介）.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 13.3 Introduction to statsmodels（statsmodels简介）\n",
    "\n",
    "[statsmodels](http://www.statsmodels.org/stable/index.html)是一个有很多统计模型的python库，能完成很多统计测试，数据探索以及可视化。它也包含一些经典的统计方法，比如贝叶斯方法和一个机器学习的模型。\n",
    "\n",
    "statsmodels中的模型包括：\n",
    "\n",
    "- 线性模型（linear models），广义线性模型（generalized linear models），鲁棒线性模型（robust linear models）\n",
    "- 线性混合效应模型（Linear mixed effects models）\n",
    "- 方差分析(ANOVA)方法（Analysis of variance (ANOVA) methods）\n",
    "- 时间序列处理（Time series processes）和状态空间模型（state space models）\n",
    "- 广义矩估计方法（Generalized method of moments）\n",
    "\n",
    "接下来我们用一些statsmodels中的工具，并了解如何使用Patsy公式和pandas DataFrame进行建模。\n",
    "\n",
    "# 1 Estimating Linear Models（估计线性模型）\n",
    "\n",
    "statsmodels中的线性模型大致分为两种：基于数组的（array-based），和基于公式的（formula-based）。调用的模块为：\n",
    "\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import statsmodels.api as sm \n",
    "import statsmodels.formula.api as smf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了演示如何使用，我们对一些随机数据生成一个线性模型："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def dnorm(mean, variance, size=1):\n",
    "    if isinstance(size, int):\n",
    "        size = size\n",
    "    return mean + np.sqrt(variance) * np.random.randn(size)\n",
    "    \n",
    "# For reproducibility \n",
    "np.random.seed(12345)\n",
    "\n",
    "N = 100 \n",
    "X = np.c_[dnorm(0, 0.4, size=N), \n",
    "          dnorm(0, 0.6, size=N), \n",
    "          dnorm(0, 0.2, size=N)] \n",
    "\n",
    "eps = dnorm(0, 0.1, size=N) \n",
    "beta = [0.1, 0.3, 0.5]\n",
    "\n",
    "y = np.dot(X, beta) + eps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(100, 3)\n",
      "(100,)\n"
     ]
    }
   ],
   "source": [
    "print(X.shape)\n",
    "print(eps.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "真正的模型用的参数是beta，dnorm的功能是产生指定平均值和方差的随机离散数据，得到："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.12946849, -1.21275292,  0.50422488],\n",
       "       [ 0.30291036, -0.43574176, -0.25417986],\n",
       "       [-0.32852189, -0.02530153,  0.13835097],\n",
       "       [-0.35147471, -0.71960511, -0.25821463],\n",
       "       [ 1.2432688 , -0.37379916, -0.52262905]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.42786349, -0.67348041, -0.09087764, -0.48949442, -0.12894109])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一个线性模型通常会有一个截距，这里我们用sm.add_constant函数添加一个截距列给X："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        , -0.12946849, -1.21275292,  0.50422488],\n",
       "       [ 1.        ,  0.30291036, -0.43574176, -0.25417986],\n",
       "       [ 1.        , -0.32852189, -0.02530153,  0.13835097],\n",
       "       [ 1.        , -0.35147471, -0.71960511, -0.25821463],\n",
       "       [ 1.        ,  1.2432688 , -0.37379916, -0.52262905]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_model = sm.add_constant(X)\n",
    "X_model[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "sm.OLS可以拟合（fit）普通最小二乘线性回归："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model = sm.OLS(y, X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "fit方法返回的是一个回顾结果对象，包含预测模型的参数和其他一些诊断数据："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.17826108,  0.22303962,  0.50095093])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results = model.fit()\n",
    "results.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在results上调用summary方法，可能得到一些详细的诊断数据："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "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.430</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.413</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   24.42</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Mon, 11 Dec 2017</td> <th>  Prob (F-statistic):</th> <td>7.44e-12</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>00:01:30</td>     <th>  Log-Likelihood:    </th> <td> -34.305</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   100</td>      <th>  AIC:               </th> <td>   74.61</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    97</td>      <th>  BIC:               </th> <td>   82.42</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     3</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>x1</th> <td>    0.1783</td> <td>    0.053</td> <td>    3.364</td> <td> 0.001</td> <td>    0.073</td> <td>    0.283</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x2</th> <td>    0.2230</td> <td>    0.046</td> <td>    4.818</td> <td> 0.000</td> <td>    0.131</td> <td>    0.315</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x3</th> <td>    0.5010</td> <td>    0.080</td> <td>    6.237</td> <td> 0.000</td> <td>    0.342</td> <td>    0.660</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 4.662</td> <th>  Durbin-Watson:     </th> <td>   2.201</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.097</td> <th>  Jarque-Bera (JB):  </th> <td>   4.098</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.481</td> <th>  Prob(JB):          </th> <td>   0.129</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 3.243</td> <th>  Cond. No.          </th> <td>    1.74</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   R-squared:                       0.430\n",
       "Model:                            OLS   Adj. R-squared:                  0.413\n",
       "Method:                 Least Squares   F-statistic:                     24.42\n",
       "Date:                Mon, 11 Dec 2017   Prob (F-statistic):           7.44e-12\n",
       "Time:                        00:01:30   Log-Likelihood:                -34.305\n",
       "No. Observations:                 100   AIC:                             74.61\n",
       "Df Residuals:                      97   BIC:                             82.42\n",
       "Df Model:                           3                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "x1             0.1783      0.053      3.364      0.001       0.073       0.283\n",
       "x2             0.2230      0.046      4.818      0.000       0.131       0.315\n",
       "x3             0.5010      0.080      6.237      0.000       0.342       0.660\n",
       "==============================================================================\n",
       "Omnibus:                        4.662   Durbin-Watson:                   2.201\n",
       "Prob(Omnibus):                  0.097   Jarque-Bera (JB):                4.098\n",
       "Skew:                           0.481   Prob(JB):                        0.129\n",
       "Kurtosis:                       3.243   Cond. No.                         1.74\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "参数的名字通常为x1, x2，以此类推。假设所有的模型参数都在一个DataFrame里："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = pd.DataFrame(X, columns=['col0', 'col1', 'col2'])\n",
    "data['y'] = y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>col0</th>\n",
       "      <th>col1</th>\n",
       "      <th>col2</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-0.129468</td>\n",
       "      <td>-1.212753</td>\n",
       "      <td>0.504225</td>\n",
       "      <td>0.427863</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.302910</td>\n",
       "      <td>-0.435742</td>\n",
       "      <td>-0.254180</td>\n",
       "      <td>-0.673480</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-0.328522</td>\n",
       "      <td>-0.025302</td>\n",
       "      <td>0.138351</td>\n",
       "      <td>-0.090878</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.351475</td>\n",
       "      <td>-0.719605</td>\n",
       "      <td>-0.258215</td>\n",
       "      <td>-0.489494</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.243269</td>\n",
       "      <td>-0.373799</td>\n",
       "      <td>-0.522629</td>\n",
       "      <td>-0.128941</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       col0      col1      col2         y\n",
       "0 -0.129468 -1.212753  0.504225  0.427863\n",
       "1  0.302910 -0.435742 -0.254180 -0.673480\n",
       "2 -0.328522 -0.025302  0.138351 -0.090878\n",
       "3 -0.351475 -0.719605 -0.258215 -0.489494\n",
       "4  1.243269 -0.373799 -0.522629 -0.128941"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在我们可以使用statsmodels formula API（公式API）和Patsy的公式字符串："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept    0.033559\n",
       "col0         0.176149\n",
       "col1         0.224826\n",
       "col2         0.514808\n",
       "dtype: float64"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results = smf.ols('y ~ col0 + col1 + col2', data=data).fit()\n",
    "results.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept    0.952188\n",
       "col0         3.319754\n",
       "col1         4.850730\n",
       "col2         6.303971\n",
       "dtype: float64"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.tvalues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以看到statsmodel返回的结果是Series，而Series的索引部分是DataFrame的列名。当我们使用公式和pandas对象的时候，不需要使用add_constant。\n",
    "\n",
    "如果得到新的数据，我们可以用预测模型的参数来进行预测："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0   -0.002327\n",
       "1   -0.141904\n",
       "2    0.041226\n",
       "3   -0.323070\n",
       "4   -0.100535\n",
       "dtype: float64"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.predict(data[:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其他一些分析、诊断、可视化工具可以自己多尝试。\n",
    "\n",
    "# 2 Estimating Time Series Processes（预测时序过程）\n",
    "\n",
    "statsmodels中的另一个类是用于时间序列分析的，其中有自动回归处理（autoregressive processes）， 卡尔曼滤波（Kalman filtering），状态空间模型（state space models），多元回归模型（multivariate autoregressive models）。\n",
    "\n",
    "让我们用自动回归结果和噪音模拟一个时间序列数据：\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "init_x = 4\n",
    "\n",
    "import random\n",
    "values = [init_x, init_x]\n",
    "N = 1000\n",
    "\n",
    "b0 = 0.8\n",
    "b1 = -0.4\n",
    "noise = dnorm(0, 0.1, N)\n",
    "for i in range(N):\n",
    "    new_x = values[-1] * b0 + values[-2] * b1 + noise[i]\n",
    "    values.append(new_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[4,\n",
       " 4,\n",
       " 1.8977509636904242,\n",
       " 0.086865262206104243,\n",
       " -0.57694691325353353,\n",
       " -0.49950238023089472]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "values[:6]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这种数据有AR(2)结构（two lags，延迟两期），延迟参数是0.8和-0.4。当我们拟合一个AR模型，我们可能不知道延迟的期间是多少，所以可以在拟合时设一个比较大的延迟数字："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "MAXLAGS = 5\n",
    "model = sm.tsa.AR(values)\n",
    "results = model.fit(MAXLAGS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "结果里的预测参数，第一个是解决，之后是两个延迟（lags）："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.00616093,  0.78446347, -0.40847891, -0.01364148,  0.01496872,\n",
       "        0.01429462])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "关于模型的更多细节，可以查看文档。"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [py35]",
   "language": "python",
   "name": "Python [py35]"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}