13.3 statsmodels介绍

statsmodels是Python进行拟合多种统计模型、进行统计试验和数据探索可视化的库。Statsmodels包含许多经典的统计方法，但没有贝叶斯方法和机器学习模型。

statsmodels包含的模型有：

  • 线性模型，广义线性模型和健壮线性模型
  • 线性混合效应模型
  • 方差（ANOVA）方法分析
  • 时间序列过程和状态空间模型
  • 广义矩估计

下面，我会使用一些基本的statsmodels工具，探索Patsy公式和pandasDataFrame对象如何使用模型接口。

估计线性模型

statsmodels有多种线性回归模型，包括从基本（比如普通最小二乘）到复杂（比如迭代加权最小二乘法）的。

statsmodels的线性模型有两种不同的接口：基于数组和基于公式。它们可以通过API模块引入：

  1. import statsmodels.api as sm
  2. import statsmodels.formula.api as smf

为了展示它们的使用方法，我们从一些随机数据生成一个线性模型：

  1. def dnorm(mean, variance, size=1):
  2.     if isinstance(size, int):
  3.         size = size,
  4.     return mean + np.sqrt(variance) * np.random.randn(*size)
  6. # For reproducibility
  7. np.random.seed(12345)
  9. N = 100
  10. X = np.c_[dnorm(0, 0.4, size=N),
  11.           dnorm(0, 0.6, size=N),
  12.           dnorm(0, 0.2, size=N)]
  13. eps = dnorm(0, 0.1, size=N)
  14. beta = [0.1, 0.3, 0.5]
  16. y = np.dot(X, beta) + eps

这里，我使用了“真实”模型和可知参数beta。此时，dnorm可用来生成正态分布数据，带有特定均值和方差。现在有：

  1. In [66]: X[:5]
  2. Out[66]: 
  3. array([[-0.1295, -1.2128,  0.5042],
  4.        [ 0.3029, -0.4357, -0.2542],
  5.        [-0.3285, -0.0253,  0.1384],
  6.        [-0.3515, -0.7196, -0.2582],
  7.        [ 1.2433, -0.3738, -0.5226]])
  9. In [67]: y[:5]
  10. Out[67]: array([ 0.4279, -0.6735, -0.0909, -0.4895,-0.1289])

像之前Patsy看到的，线性模型通常要拟合一个截距。sm.add_constant函数可以添加一个截距的列到现存的矩阵：

  1. In [68]: X_model = sm.add_constant(X)
  3. In [69]: X_model[:5]
  4. Out[69]: 
  5. array([[ 1.    , -0.1295, -1.2128,  0.5042],
  6.        [ 1.    ,  0.3029, -0.4357, -0.2542],
  7.        [ 1.    , -0.3285, -0.0253,  0.1384],
  8.        [ 1.    , -0.3515, -0.7196, -0.2582],
  9.        [ 1.    ,  1.2433, -0.3738, -0.5226]])

sm.OLS类可以拟合一个普通最小二乘回归：

  1. In [70]: model = sm.OLS(y, X)

这个模型的fit方法返回了一个回归结果对象，它包含估计的模型参数和其它内容：

  1. In [71]: results = model.fit()
  3. In [72]: results.params
  4. Out[72]: array([ 0.1783,  0.223 ,  0.501 ])

对结果使用summary方法可以打印模型的详细诊断结果：

  1. In [73]: print(results.summary())
  2. OLS Regression Results                            
  3. ==============================================================================
  4. Dep. Variable:                      y   R-squared:                       0.430
  5. Model:                            OLS   Adj. R-squared:                  0.413
  6. Method:                 Least Squares   F-statistic:                     24.42
  7. Date:                Mon, 25 Sep 2017   Prob (F-statistic):           7.44e-12
  8. Time:                        14:06:15   Log-Likelihood:                -34.305
  9. No. Observations:                 100   AIC:                             74.61
  10. Df Residuals:                      97   BIC:                             82.42
  11. Df Model:                           3                                         
  12. Covariance Type:            nonrobust                                         
  13. ==============================================================================
  14.                  coef    std err          t      P>|t|      [0.025      0.975]
  15. ------------------------------------------------------------------------------
  16. x1             0.1783      0.053      3.364      0.001       0.073       0.283
  17. x2             0.2230      0.046      4.818      0.000       0.131       0.315
  18. x3             0.5010      0.080      6.237      0.000       0.342       0.660
  19. ==============================================================================
  20. Omnibus:                        4.662   Durbin-Watson:                   2.201
  21. Prob(Omnibus):                  0.097   Jarque-Bera (JB):                4.098
  22. Skew:                           0.481   Prob(JB):                        0.129
  23. Kurtosis:                       3.243   Cond. No.
  24. 1.74
  25. ==============================================================================
  26. Warnings:
  27. [1] Standard Errors assume that the covariance matrix of the errors is correctly 
  28. specified.

这里的参数名为通用名x1, x2等等。假设所有的模型参数都在一个DataFrame中：

  1. In [74]: data = pd.DataFrame(X, columns=['col0', 'col1', 'col2'])
  3. In [75]: data['y'] = y
  5. In [76]: data[:5]
  6. Out[76]: 
  7.        col0      col1      col2         y
  8. 0 -0.129468 -1.212753  0.504225  0.427863
  9. 1  0.302910 -0.435742 -0.254180 -0.673480
  10. 2 -0.328522 -0.025302  0.138351 -0.090878
  11. 3 -0.351475 -0.719605 -0.258215 -0.489494
  12. 4  1.243269 -0.373799 -0.522629 -0.128941

现在，我们使用statsmodels的公式API和Patsy的公式字符串：

  1. In [77]: results = smf.ols('y ~ col0 + col1 + col2', data=data).fit()
  3. In [78]: results.params
  4. Out[78]: 
  5. Intercept    0.033559
  6. col0         0.176149
  7. col1         0.224826
  8. col2         0.514808
  9. dtype: float64
  11. In [79]: results.tvalues
  12. Out[79]: 
  13. Intercept    0.952188
  14. col0         3.319754
  15. col1         4.850730
  16. col2         6.303971
  17. dtype: float64

观察下statsmodels是如何返回Series结果的，附带有DataFrame的列名。当使用公式和pandas对象时，我们不需要使用add_constant。

给出一个样本外数据，你可以根据估计的模型参数计算预测值：

  1. In [80]: results.predict(data[:5])
  2. Out[80]: 
  3. 0   -0.002327
  4. 1   -0.141904
  5. 2    0.041226
  6. 3   -0.323070
  7. 4   -0.100535
  8. dtype: float64

statsmodels的线性模型结果还有其它的分析、诊断和可视化工具。除了普通最小二乘模型，还有其它的线性模型。

估计时间序列过程

statsmodels的另一模型类是进行时间序列分析，包括自回归过程、卡尔曼滤波和其它态空间模型，和多元自回归模型。

用自回归结构和噪声来模拟一些时间序列数据：

  1. init_x = 4
  3. import random
  4. values = [init_x, init_x]
  5. N = 1000
  7. b0 = 0.8
  8. b1 = -0.4
  9. noise = dnorm(0, 0.1, N)
  10. for i in range(N):
  11.     new_x = values[-1] * b0 + values[-2] * b1 + noise[i]
  12.     values.append(new_x)

这个数据有AR(2)结构（两个延迟），参数是0.8和-0.4。拟合AR模型时，你可能不知道滞后项的个数，因此可以用较多的滞后量来拟合这个模型：

  1. In [82]: MAXLAGS = 5
  3. In [83]: model = sm.tsa.AR(values)
  5. In [84]: results = model.fit(MAXLAGS)

结果中的估计参数首先是截距，其次是前两个参数的估计值：

  1. In [85]: results.params
  2. Out[85]: array([-0.0062,  0.7845, -0.4085, -0.0136,  0.015 ,  0.0143])

更多的细节以及如何解释结果超出了本书的范围，可以通过statsmodels文档学习更多。