Python 数据分析（PYDA）第三版（六）（1）https://developer.aliyun.com/article/1482396
12.3 statsmodels 简介
statsmodels是一个用于拟合许多种统计模型、执行统计检验以及数据探索和可视化的 Python 库。statsmodels 包含更多“经典”的频率统计方法，而贝叶斯方法和机器学习模型则在其他库中找到。
在 statsmodels 中找到的一些模型类型包括：
- 线性模型、广义线性模型和鲁棒线性模型
- 线性混合效应模型
- 方差分析（ANOVA）方法
- 时间序列过程和状态空间模型
- 广义矩估计法
在接下来的几页中，我们将使用 statsmodels 中的一些基本工具，并探索如何使用 Patsy 公式和 pandas DataFrame 对象的建模接口。如果您之前在 Patsy 讨论中没有安装 statsmodels，现在可以使用以下命令进行安装：
conda install statsmodels
估计线性模型
statsmodels 中有几种线性回归模型，从更基本的（例如普通最小二乘法）到更复杂的（例如迭代重新加权最小二乘法）。
statsmodels 中的线性模型有两种不同的主要接口：基于数组和基于公式。可以通过以下 API 模块导入来访问这些接口：
import statsmodels.api as sm import statsmodels.formula.api as smf
为了展示如何使用这些方法，我们从一些随机数据生成一个线性模型。在 Jupyter 中运行以下代码：
# To make the example reproducible rng = np.random.default_rng(seed=12345) def dnorm(mean, variance, size=1): if isinstance(size, int): size = size, return mean + np.sqrt(variance) * rng.standard_normal(*size) N = 100 X = np.c_[dnorm(0, 0.4, size=N), dnorm(0, 0.6, size=N), dnorm(0, 0.2, size=N)] eps = dnorm(0, 0.1, size=N) beta = [0.1, 0.3, 0.5] y = np.dot(X, beta) + eps
在这里，我写下了具有已知参数beta的“真实”模型。在这种情况下，dnorm是一个用于生成具有特定均值和方差的正态分布数据的辅助函数。现在我们有：
In [66]: X[:5] Out[66]: array([[-0.9005, -0.1894, -1.0279], [ 0.7993, -1.546 , -0.3274], [-0.5507, -0.1203, 0.3294], [-0.1639, 0.824 , 0.2083], [-0.0477, -0.2131, -0.0482]]) In [67]: y[:5] Out[67]: array([-0.5995, -0.5885, 0.1856, -0.0075, -0.0154])
通常使用截距项拟合线性模型，就像我们之前在 Patsy 中看到的那样。sm.add_constant函数可以向现有矩阵添加一个截距列：
In [68]: X_model = sm.add_constant(X) In [69]: X_model[:5] Out[69]: array([[ 1. , -0.9005, -0.1894, -1.0279], [ 1. , 0.7993, -1.546 , -0.3274], [ 1. , -0.5507, -0.1203, 0.3294], [ 1. , -0.1639, 0.824 , 0.2083], [ 1. , -0.0477, -0.2131, -0.0482]])
sm.OLS类可以拟合普通最小二乘线性回归：
In [70]: model = sm.OLS(y, X)
模型的fit方法返回一个包含估计模型参数和其他诊断信息的回归结果对象：
In [71]: results = model.fit() In [72]: results.params Out[72]: array([0.0668, 0.268 , 0.4505])
results上的summary方法可以打印出模型的诊断输出：
In [73]: print(results.summary()) OLS Regression Results ================================================================================= ====== Dep. Variable: y R-squared (uncentered): 0.469 Model: OLS Adj. R-squared (uncentered): 0.452 Method: Least Squares F-statistic: 28.51 Date: Wed, 12 Apr 2023 Prob (F-statistic): 2. 66e-13 Time: 13:09:20 Log-Likelihood: - 25.611 No. Observations: 100 AIC: 57.22 Df Residuals: 97 BIC: 65.04 Df Model: 3 Covariance Type: nonrobust ============================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------ x1 0.0668 0.054 1.243 0.217 -0.040 0.174 x2 0.2680 0.042 6.313 0.000 0.184 0.352 x3 0.4505 0.068 6.605 0.000 0.315 0.586 ============================================================================== Omnibus: 0.435 Durbin-Watson: 1.869 Prob(Omnibus): 0.805 Jarque-Bera (JB): 0.301 Skew: 0.134 Prob(JB): 0.860 Kurtosis: 2.995 Cond. No. 1.64 ============================================================================== Notes: [1] R² is computed without centering (uncentered) since the model does not contai n a constant. [2] Standard Errors assume that the covariance matrix of the errors is correctly specified.
这里的参数名称已经被赋予了通用名称x1, x2等。假设所有模型参数都在一个 DataFrame 中：
In [74]: data = pd.DataFrame(X, columns=['col0', 'col1', 'col2']) In [75]: data['y'] = y In [76]: data[:5] Out[76]: col0 col1 col2 y 0 -0.900506 -0.189430 -1.027870 -0.599527 1 0.799252 -1.545984 -0.327397 -0.588454 2 -0.550655 -0.120254 0.329359 0.185634 3 -0.163916 0.824040 0.208275 -0.007477 4 -0.047651 -0.213147 -0.048244 -0.015374
现在我们可以使用 statsmodels 的公式 API 和 Patsy 公式字符串：
In [77]: results = smf.ols('y ~ col0 + col1 + col2', data=data).fit() In [78]: results.params Out[78]: Intercept -0.020799 col0 0.065813 col1 0.268970 col2 0.449419 dtype: float64 In [79]: results.tvalues Out[79]: Intercept -0.652501 col0 1.219768 col1 6.312369 col2 6.567428 dtype: float64
注意 statsmodels 如何将结果返回为带有 DataFrame 列名称附加的 Series。在使用公式和 pandas 对象时，我们也不需要使用add_constant。
给定新的样本外数据，可以根据估计的模型参数计算预测值：
In [80]: results.predict(data[:5]) Out[80]: 0 -0.592959 1 -0.531160 2 0.058636 3 0.283658 4 -0.102947 dtype: float64
在 statsmodels 中有许多用于分析、诊断和可视化线性模型结果的附加工具，您可以探索。除了普通最小二乘法之外，还有其他类型的线性模型。
估计时间序列过程
statsmodels 中的另一类模型是用于时间序列分析的模型。其中包括自回归过程、卡尔曼滤波和其他状态空间模型以及多变量自回归模型。
让我们模拟一些具有自回归结构和噪声的时间序列数据。在 Jupyter 中运行以下代码：
init_x = 4 values = [init_x, init_x] N = 1000 b0 = 0.8 b1 = -0.4 noise = dnorm(0, 0.1, N) for i in range(N): new_x = values[-1] * b0 + values[-2] * b1 + noise[i] values.append(new_x)
这个数据具有 AR(2)结构（两个滞后），参数为0.8和-0.4。当拟合 AR 模型时，您可能不知道要包括的滞后项的数量，因此可以使用一些更大数量的滞后项来拟合模型：
In [82]: from statsmodels.tsa.ar_model import AutoReg In [83]: MAXLAGS = 5 In [84]: model = AutoReg(values, MAXLAGS) In [85]: results = model.fit()
结果中的估计参数首先是截距，接下来是前两个滞后的估计值：
In [86]: results.params Out[86]: array([ 0.0235, 0.8097, -0.4287, -0.0334, 0.0427, -0.0567])
这些模型的更深层细节以及如何解释它们的结果超出了我在本书中可以涵盖的范围，但在 statsmodels 文档中还有很多内容等待探索。
12.4 scikit-learn 简介
scikit-learn是最广泛使用和信任的通用 Python 机器学习工具包之一。它包含广泛的标准监督和无监督机器学习方法，具有模型选择和评估工具，数据转换，数据加载和模型持久性。这些模型可用于分类，聚类，预测和其他常见任务。您可以像这样从 conda 安装 scikit-learn：
conda install scikit-learn
有很多在线和印刷资源可供学习机器学习以及如何应用类似 scikit-learn 的库来解决实际问题。在本节中，我将简要介绍 scikit-learn API 风格。
scikit-learn 中的 pandas 集成在近年来显著改善，当您阅读本文时，它可能已经进一步改进。我鼓励您查看最新的项目文档。
作为本章的示例，我使用了一份来自 Kaggle 竞赛的经典数据集，关于 1912 年泰坦尼克号上乘客生存率。我们使用 pandas 加载训练和测试数据集：
In [87]: train = pd.read_csv('datasets/titanic/train.csv') In [88]: test = pd.read_csv('datasets/titanic/test.csv') In [89]: train.head(4) Out[89]: PassengerId Survived Pclass 0 1 0 3 \ 1 2 1 1 2 3 1 3 3 4 1 1 Name Sex Age SibSp 0 Braund, Mr. Owen Harris male 22.0 1 \ 1 Cumings, Mrs. John Bradley (Florence Briggs Thayer) female 38.0 1 2 Heikkinen, Miss. Laina female 26.0 0 3 Futrelle, Mrs. Jacques Heath (Lily May Peel) female 35.0 1 Parch Ticket Fare Cabin Embarked 0 0 A/5 21171 7.2500 NaN S 1 0 PC 17599 71.2833 C85 C 2 0 STON/O2\. 3101282 7.9250 NaN S 3 0 113803 53.1000 C123 S
像 statsmodels 和 scikit-learn 这样的库通常无法处理缺失数据，因此我们查看列，看看是否有包含缺失数据的列：
In [90]: train.isna().sum() Out[90]: PassengerId 0 Survived 0 Pclass 0 Name 0 Sex 0 Age 177 SibSp 0 Parch 0 Ticket 0 Fare 0 Cabin 687 Embarked 2 dtype: int64 In [91]: test.isna().sum() Out[91]: PassengerId 0 Pclass 0 Name 0 Sex 0 Age 86 SibSp 0 Parch 0 Ticket 0 Fare 1 Cabin 327 Embarked 0 dtype: int64
在统计学和机器学习的示例中，一个典型的任务是根据数据中的特征预测乘客是否会生存。模型在训练数据集上拟合，然后在外样本测试数据集上进行评估。
我想使用Age作为预测变量，但它有缺失数据。有很多方法可以进行缺失数据插补，但我将使用训练数据集的中位数来填充两个表中的空值：
In [92]: impute_value = train['Age'].median() In [93]: train['Age'] = train['Age'].fillna(impute_value) In [94]: test['Age'] = test['Age'].fillna(impute_value)
现在我们需要指定我们的模型。我添加一个名为IsFemale的列，作为'Sex'列的编码版本：
In [95]: train['IsFemale'] = (train['Sex'] == 'female').astype(int) In [96]: test['IsFemale'] = (test['Sex'] == 'female').astype(int)
然后我们决定一些模型变量并创建 NumPy 数组：
In [97]: predictors = ['Pclass', 'IsFemale', 'Age'] In [98]: X_train = train[predictors].to_numpy() In [99]: X_test = test[predictors].to_numpy() In [100]: y_train = train['Survived'].to_numpy() In [101]: X_train[:5] Out[101]: array([[ 3., 0., 22.], [ 1., 1., 38.], [ 3., 1., 26.], [ 1., 1., 35.], [ 3., 0., 35.]]) In [102]: y_train[:5] Out[102]: array([0, 1, 1, 1, 0])
我不断言这是一个好模型或这些特征是否被正确设计。我们使用 scikit-learn 中的LogisticRegression模型并创建一个模型实例：
In [103]: from sklearn.linear_model import LogisticRegression In [104]: model = LogisticRegression()
我们可以使用模型的fit方法将此模型拟合到训练数据中：
In [105]: model.fit(X_train, y_train) Out[105]: LogisticRegression()
现在，我们可以使用model.predict为测试数据集进行预测：
In [106]: y_predict = model.predict(X_test) In [107]: y_predict[:10] Out[107]: array([0, 0, 0, 0, 1, 0, 1, 0, 1, 0])
如果您有测试数据集的真实值，可以计算准确率百分比或其他错误度量：
(y_true == y_predict).mean()
在实践中，模型训练通常存在许多额外的复杂层。许多模型具有可以调整的参数，并且有一些技术，如交叉验证可用于参数调整，以避免过度拟合训练数据。这通常可以提供更好的预测性能或对新数据的鲁棒性。
交叉验证通过拆分训练数据来模拟外样本预测。根据像均方误差这样的模型准确度得分，您可以对模型参数执行网格搜索。一些模型，如逻辑回归，具有内置交叉验证的估计器类。例如，LogisticRegressionCV类可以与一个参数一起使用，该参数指示在模型正则化参数C上执行多精细的网格搜索：
In [108]: from sklearn.linear_model import LogisticRegressionCV In [109]: model_cv = LogisticRegressionCV(Cs=10) In [110]: model_cv.fit(X_train, y_train) Out[110]: LogisticRegressionCV()
手动进行交叉验证，可以使用cross_val_score辅助函数，该函数处理数据拆分过程。例如，要对我们的模型进行四个不重叠的训练数据拆分进行交叉验证，我们可以这样做：
In [111]: from sklearn.model_selection import cross_val_score In [112]: model = LogisticRegression(C=10) In [113]: scores = cross_val_score(model, X_train, y_train, cv=4) In [114]: scores Out[114]: array([0.7758, 0.7982, 0.7758, 0.7883])
默认的评分指标取决于模型，但可以选择一个明确的评分函数。交叉验证模型训练时间较长，但通常可以获得更好的模型性能。
12.5 结论
虽然我只是浅尝了一些 Python 建模库的表面，但有越来越多的框架适用于各种统计和机器学习，要么是用 Python 实现的，要么有 Python 用户界面。
这本书专注于数据整理，但还有许多其他专门用于建模和数据科学工具的书籍。一些优秀的书籍包括：
- 《Python 机器学习入门》作者 Andreas Müller 和 Sarah Guido（O’Reilly）
- 《Python 数据科学手册》作者 Jake VanderPlas（O’Reilly）
- 《从零开始的数据科学：Python 基础》作者 Joel Grus（O’Reilly）
- 《Python 机器学习》作者 Sebastian Raschka 和 Vahid Mirjalili（Packt Publishing）
- 《使用 Scikit-Learn、Keras 和 TensorFlow 进行实践机器学习》作者 Aurélien Géron（O’Reilly）
尽管书籍可以是学习的宝贵资源，但当底层的开源软件发生变化时，它们有时会变得过时。熟悉各种统计或机器学习框架的文档是一个好主意，以便了解最新功能和 API。
Python 数据分析（PYDA）第三版（六）（3）https://developer.aliyun.com/article/1482398
