原文：
wesmckinney.com/book/译者：飞龙
协议：CC BY-NC-SA 4.0
十二、Python 建模库介绍
原文：
wesmckinney.com/book/modeling译者：飞龙
协议：CC BY-NC-SA 4.0
此开放访问网络版本的《Python 数据分析第三版》现已作为印刷版和数字版的伴侣提供。如果您发现任何勘误，请在此处报告。请注意，由 Quarto 生成的本站点的某些方面与 O’Reilly 的印刷版和电子书版本的格式不同。
如果您发现本书的在线版本有用，请考虑订购纸质版或无 DRM 的电子书以支持作者。本网站的内容不得复制或再生产。代码示例采用 MIT 许可，可在 GitHub 或 Gitee 上找到。
在本书中，我专注于为在 Python 中进行数据分析提供编程基础。由于数据分析师和科学家经常报告花费大量时间进行数据整理和准备，因此本书的结构反映了掌握这些技术的重要性。
您用于开发模型的库将取决于应用程序。许多统计问题可以通过简单的技术解决，如普通最小二乘回归，而其他问题可能需要更高级的机器学习方法。幸运的是，Python 已经成为实现分析方法的首选语言之一，因此在完成本书后，您可以探索许多工具。
在本章中，我将回顾一些 pandas 的特性，这些特性在您在 pandas 中进行数据整理和模型拟合和评分之间来回切换时可能会有所帮助。然后，我将简要介绍两个流行的建模工具包，statsmodels和scikit-learn。由于这两个项目都足够庞大，值得有自己的专门书籍，因此我没有尝试全面介绍，而是建议您查阅这两个项目的在线文档，以及一些其他基于 Python 的数据科学、统计学和机器学习书籍。
12.1 pandas 与模型代码之间的接口
模型开发的常见工作流程是使用 pandas 进行数据加载和清理，然后切换到建模库来构建模型本身。模型开发过程中的一个重要部分被称为特征工程，在机器学习中。这可以描述从原始数据集中提取信息的任何数据转换或分析，这些信息在建模环境中可能有用。我们在本书中探讨的数据聚合和 GroupBy 工具经常在特征工程环境中使用。
虽然“好”的特征工程的细节超出了本书的范围，但我将展示一些方法，使在 pandas 中进行数据操作和建模之间的切换尽可能轻松。
pandas 与其他分析库之间的接触点通常是 NumPy 数组。要将 DataFrame 转换为 NumPy 数组，请使用to_numpy方法：
In [12]: data = pd.DataFrame({
....: 'x0': [1, 2, 3, 4, 5],
....: 'x1': [0.01, -0.01, 0.25, -4.1, 0.],
....: 'y': [-1.5, 0., 3.6, 1.3, -2.]})
In [13]: data
Out[13]:
x0 x1 y
0 1 0.01 -1.5
1 2 -0.01 0.0
2 3 0.25 3.6
3 4 -4.10 1.3
4 5 0.00 -2.0
In [14]: data.columns
Out[14]: Index(['x0', 'x1', 'y'], dtype='object')
In [15]: data.to_numpy()
Out[15]:
array([[ 1. , 0.01, -1.5 ],
[ 2. , -0.01, 0. ],
[ 3. , 0.25, 3.6 ],
[ 4. , -4.1 , 1.3 ],
[ 5. , 0. , -2. ]])
回到 DataFrame，正如您可能从前几章中记得的那样，您可以传递一个二维的 ndarray，其中包含可选的列名：
In [16]: df2 = pd.DataFrame(data.to_numpy(), columns=['one', 'two', 'three'])
In [17]: df2
Out[17]:
one two three
0 1.0 0.01 -1.5
1 2.0 -0.01 0.0
2 3.0 0.25 3.6
3 4.0 -4.10 1.3
4 5.0 0.00 -2.0
to_numpy方法旨在在数据是同质的情况下使用，例如所有的数值类型。如果您有异构数据，结果将是一个 Python 对象的 ndarray：
In [18]: df3 = data.copy()
In [19]: df3['strings'] = ['a', 'b', 'c', 'd', 'e']
In [20]: df3
Out[20]:
x0 x1 y strings
0 1 0.01 -1.5 a
1 2 -0.01 0.0 b
2 3 0.25 3.6 c
3 4 -4.10 1.3 d
4 5 0.00 -2.0 e
In [21]: df3.to_numpy()
Out[21]:
array([[1, 0.01, -1.5, 'a'],
[2, -0.01, 0.0, 'b'],
[3, 0.25, 3.6, 'c'],
[4, -4.1, 1.3, 'd'],
[5, 0.0, -2.0, 'e']], dtype=object)
对于某些模型，您可能希望仅使用部分列。我建议使用loc索引和to_numpy：
In [22]: model_cols = ['x0', 'x1']
In [23]: data.loc[:, model_cols].to_numpy()
Out[23]:
array([[ 1. , 0.01],
[ 2. , -0.01],
[ 3. , 0.25],
[ 4. , -4.1 ],
[ 5. , 0. ]])
一些库原生支持 pandas，并自动完成一些工作：从 DataFrame 转换为 NumPy，并将模型参数名称附加到输出表或 Series 的列上。在其他情况下，您将不得不手动执行这种“元数据管理”。
在 Ch 7.5：分类数据中，我们看过 pandas 的Categorical类型和pandas.get_dummies函数。假设我们的示例数据集中有一个非数字列：
In [24]: data['category'] = pd.Categorical(['a', 'b', 'a', 'a', 'b'],
....: categories=['a', 'b'])
In [25]: data
Out[25]:
x0 x1 y category
0 1 0.01 -1.5 a
1 2 -0.01 0.0 b
2 3 0.25 3.6 a
3 4 -4.10 1.3 a
4 5 0.00 -2.0 b
如果我们想用虚拟变量替换'category'列，我们创建虚拟变量，删除'category'列，然后将结果连接：
In [26]: dummies = pd.get_dummies(data.category, prefix='category',
....: dtype=float)
In [27]: data_with_dummies = data.drop('category', axis=1).join(dummies)
In [28]: data_with_dummies
Out[28]:
x0 x1 y category_a category_b
0 1 0.01 -1.5 1.0 0.0
1 2 -0.01 0.0 0.0 1.0
2 3 0.25 3.6 1.0 0.0
3 4 -4.10 1.3 1.0 0.0
4 5 0.00 -2.0 0.0 1.0
使用虚拟变量拟合某些统计模型时存在一些微妙之处。当您拥有不仅仅是简单数字列时，使用 Patsy（下一节的主题）可能更简单且更不容易出错。
12.2 使用 Patsy 创建模型描述
Patsy是一个用于描述统计模型（尤其是线性模型）的 Python 库，它使用基于字符串的“公式语法”，受到 R 和 S 统计编程语言使用的公式语法的启发（但并非完全相同）。在安装 statsmodels 时会自动安装它：
conda install statsmodels
Patsy 在为 statsmodels 指定线性模型方面得到很好的支持，因此我将重点介绍一些主要功能，以帮助您快速上手。Patsy 的公式是一种特殊的字符串语法，看起来像：
y ~ x0 + x1
语法a + b并不意味着将a加到b，而是这些是为模型创建的设计矩阵中的项。patsy.dmatrices函数接受一个公式字符串以及一个数据集（可以是 DataFrame 或数组字典），并为线性模型生成设计矩阵：
In [29]: data = pd.DataFrame({
....: 'x0': [1, 2, 3, 4, 5],
....: 'x1': [0.01, -0.01, 0.25, -4.1, 0.],
....: 'y': [-1.5, 0., 3.6, 1.3, -2.]})
In [30]: data
Out[30]:
x0 x1 y
0 1 0.01 -1.5
1 2 -0.01 0.0
2 3 0.25 3.6
3 4 -4.10 1.3
4 5 0.00 -2.0
In [31]: import patsy
In [32]: y, X = patsy.dmatrices('y ~ x0 + x1', data)
现在我们有：
In [33]: y
Out[33]:
DesignMatrix with shape (5, 1)
y
-1.5
0.0
3.6
1.3
-2.0
Terms:
'y' (column 0)
In [34]: X
Out[34]:
DesignMatrix with shape (5, 3)
Intercept x0 x1
1 1 0.01
1 2 -0.01
1 3 0.25
1 4 -4.10
1 5 0.00
Terms:
'Intercept' (column 0)
'x0' (column 1)
'x1' (column 2)
这些 Patsy DesignMatrix实例是带有附加元数据的 NumPy ndarrays：
In [35]: np.asarray(y)
Out[35]:
array([[-1.5],
[ 0. ],
[ 3.6],
[ 1.3],
[-2. ]])
In [36]: np.asarray(X)
Out[36]:
array([[ 1. , 1. , 0.01],
[ 1. , 2. , -0.01],
[ 1. , 3. , 0.25],
[ 1. , 4. , -4.1 ],
[ 1. , 5. , 0. ]])
您可能会想知道Intercept项是从哪里来的。这是线性模型（如普通最小二乘回归）的一个约定。您可以通过在模型中添加+ 0项来抑制截距：
In [37]: patsy.dmatrices('y ~ x0 + x1 + 0', data)[1]
Out[37]:
DesignMatrix with shape (5, 2)
x0 x1
1 0.01
2 -0.01
3 0.25
4 -4.10
5 0.00
Terms:
'x0' (column 0)
'x1' (column 1)
Patsy 对象可以直接传递到像numpy.linalg.lstsq这样的算法中，该算法执行普通最小二乘回归：
In [38]: coef, resid, _, _ = np.linalg.lstsq(X, y, rcond=None)
模型元数据保留在design_info属性中，因此您可以重新附加模型列名称到拟合系数以获得一个 Series，例如：
In [39]: coef
Out[39]:
array([[ 0.3129],
[-0.0791],
[-0.2655]])
In [40]: coef = pd.Series(coef.squeeze(), index=X.design_info.column_names)
In [41]: coef
Out[41]:
Intercept 0.312910
x0 -0.079106
x1 -0.265464
dtype: float64
Patsy 公式中的数据转换
您可以将 Python 代码混合到您的 Patsy 公式中；在评估公式时，库将尝试在封闭范围中找到您使用的函数：
In [42]: y, X = patsy.dmatrices('y ~ x0 + np.log(np.abs(x1) + 1)', data)
In [43]: X
Out[43]:
DesignMatrix with shape (5, 3)
Intercept x0 np.log(np.abs(x1) + 1)
1 1 0.00995
1 2 0.00995
1 3 0.22314
1 4 1.62924
1 5 0.00000
Terms:
'Intercept' (column 0)
'x0' (column 1)
'np.log(np.abs(x1) + 1)' (column 2)
一些常用的变量转换包括标准化（均值为 0，方差为 1）和中心化（减去均值）。Patsy 具有内置函数用于此目的：
In [44]: y, X = patsy.dmatrices('y ~ standardize(x0) + center(x1)', data)
In [45]: X
Out[45]:
DesignMatrix with shape (5, 3)
Intercept standardize(x0) center(x1)
1 -1.41421 0.78
1 -0.70711 0.76
1 0.00000 1.02
1 0.70711 -3.33
1 1.41421 0.77
Terms:
'Intercept' (column 0)
'standardize(x0)' (column 1)
'center(x1)' (column 2)
作为建模过程的一部分，您可以在一个数据集上拟合模型，然后基于另一个数据集评估模型。这可能是一个保留部分或稍后观察到的新数据。当应用诸如中心化和标准化之类的转换时，您在使用模型基于新数据形成预测时应当小心。这些被称为有状态转换，因为在转换新数据时必须使用原始数据集的统计数据，如均值或标准差。
patsy.build_design_matrices函数可以使用原始样本内数据的保存信息对新的样本外数据应用转换：
In [46]: new_data = pd.DataFrame({
....: 'x0': [6, 7, 8, 9],
....: 'x1': [3.1, -0.5, 0, 2.3],
....: 'y': [1, 2, 3, 4]})
In [47]: new_X = patsy.build_design_matrices([X.design_info], new_data)
In [48]: new_X
Out[48]:
[DesignMatrix with shape (4, 3)
Intercept standardize(x0) center(x1)
1 2.12132 3.87
1 2.82843 0.27
1 3.53553 0.77
1 4.24264 3.07
Terms:
'Intercept' (column 0)
'standardize(x0)' (column 1)
'center(x1)' (column 2)]
因为 Patsy 公式中加号（+）并不表示加法，所以当您想按名称从数据集中添加列时，您必须将它们包装在特殊的I函数中：
In [49]: y, X = patsy.dmatrices('y ~ I(x0 + x1)', data)
In [50]: X
Out[50]:
DesignMatrix with shape (5, 2)
Intercept I(x0 + x1)
1 1.01
1 1.99
1 3.25
1 -0.10
1 5.00
Terms:
'Intercept' (column 0)
'I(x0 + x1)' (column 1)
Patsy 在patsy.builtins模块中还有几个内置转换。请查看在线文档以获取更多信息。
分类数据有一类特殊的转换，接下来我会解释。
分类数据和 Patsy
非数字数据可以以多种不同的方式转换为模型设计矩阵。本书不涉及这个主题的完整处理，最好是在统计课程中学习。
当您在 Patsy 公式中使用非数字术语时，默认情况下它们会被转换为虚拟变量。如果有一个截距，将会有一个级别被排除以避免共线性：
In [51]: data = pd.DataFrame({
....: 'key1': ['a', 'a', 'b', 'b', 'a', 'b', 'a', 'b'],
....: 'key2': [0, 1, 0, 1, 0, 1, 0, 0],
....: 'v1': [1, 2, 3, 4, 5, 6, 7, 8],
....: 'v2': [-1, 0, 2.5, -0.5, 4.0, -1.2, 0.2, -1.7]
....: })
In [52]: y, X = patsy.dmatrices('v2 ~ key1', data)
In [53]: X
Out[53]:
DesignMatrix with shape (8, 2)
Intercept key1[T.b]
1 0
1 0
1 1
1 1
1 0
1 1
1 0
1 1
Terms:
'Intercept' (column 0)
'key1' (column 1)
如果从模型中省略截距，那么每个类别值的列将包含在模型设计矩阵中：
In [54]: y, X = patsy.dmatrices('v2 ~ key1 + 0', data)
In [55]: X
Out[55]:
DesignMatrix with shape (8, 2)
key1[a] key1[b]
1 0
1 0
0 1
0 1
1 0
0 1
1 0
0 1
Terms:
'key1' (columns 0:2)
数值列可以使用C函数解释为分类列：
In [56]: y, X = patsy.dmatrices('v2 ~ C(key2)', data)
In [57]: X
Out[57]:
DesignMatrix with shape (8, 2)
Intercept C(key2)[T.1]
1 0
1 1
1 0
1 1
1 0
1 1
1 0
1 0
Terms:
'Intercept' (column 0)
'C(key2)' (column 1)
当您在模型中使用多个分类项时，情况可能会更加复杂，因为您可以包括形式为key1:key2的交互项，例如在方差分析（ANOVA）模型中使用：
In [58]: data['key2'] = data['key2'].map({0: 'zero', 1: 'one'})
In [59]: data
Out[59]:
key1 key2 v1 v2
0 a zero 1 -1.0
1 a one 2 0.0
2 b zero 3 2.5
3 b one 4 -0.5
4 a zero 5 4.0
5 b one 6 -1.2
6 a zero 7 0.2
7 b zero 8 -1.7
In [60]: y, X = patsy.dmatrices('v2 ~ key1 + key2', data)
In [61]: X
Out[61]:
DesignMatrix with shape (8, 3)
Intercept key1[T.b] key2[T.zero]
1 0 1
1 0 0
1 1 1
1 1 0
1 0 1
1 1 0
1 0 1
1 1 1
Terms:
'Intercept' (column 0)
'key1' (column 1)
'key2' (column 2)
In [62]: y, X = patsy.dmatrices('v2 ~ key1 + key2 + key1:key2', data)
In [63]: X
Out[63]:
DesignMatrix with shape (8, 4)
Intercept key1[T.b] key2[T.zero] key1[T.b]:key2[T.zero]
1 0 1 0
1 0 0 0
1 1 1 1
1 1 0 0
1 0 1 0
1 1 0 0
1 0 1 0
1 1 1 1
Terms:
'Intercept' (column 0)
'key1' (column 1)
'key2' (column 2)
'key1:key2' (column 3)
Patsy 提供了其他转换分类数据的方法，包括具有特定顺序的项的转换。有关更多信息，请参阅在线文档。
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。
十三、数据分析示例
原文：
wesmckinney.com/book/data-analysis-examples译者：飞龙
协议：CC BY-NC-SA 4.0
此开放访问网络版本的《Python 数据分析第三版》现已作为印刷版和数字版的伴侣提供。如果您发现任何勘误，请在此处报告。请注意，由 Quarto 生成的本站点的某些方面与 O’Reilly 的印刷版和电子书版本的格式不同。
如果您发现本书的在线版本有用，请考虑订购纸质版或无 DRM 的电子书以支持作者。本网站的内容不得复制或再生产。代码示例采用 MIT 许可，可在 GitHub 或 Gitee 上找到。
现在我们已经到达本书的最后一章，我们将查看一些真实世界的数据集。对于每个数据集，我们将使用本书中介绍的技术从原始数据中提取含义。演示的技术可以应用于各种其他数据集。本章包含一系列杂例数据集，您可以使用这些数据集练习本书中的工具。
示例数据集可在本书附带的GitHub 存储库中找到。如果无法访问 GitHub，还可以从Gitee 上的存储库镜像获取它们。
13.1 Bitly Data from 1.USA.gov
2011 年，URL 缩短服务Bitly与美国政府网站USA.gov合作，提供从缩短链接以*.gov或.mil*结尾的用户收集的匿名数据的源。2011 年，可下载的文本文件提供了实时数据以及每小时的快照。本文撰写时（2022 年），该服务已关闭，但我们保留了一份数据文件用于本书的示例。
在每个文件的每一行中，每小时快照包含一种称为 JSON 的常见网络数据形式，JSON 代表 JavaScript 对象表示法。例如，如果我们只读取文件的第一行，可能会看到类似于这样的内容：
In [5]: path = "datasets/bitly_usagov/example.txt"
In [6]: with open(path) as f:
...: print(f.readline())
...:
{ "a": "Mozilla\\/5.0 (Windows NT 6.1; WOW64) AppleWebKit\\/535.11
(KHTML, like Gecko) Chrome\\/17.0.963.78 Safari\\/535.11", "c": "US", "nk": 1,
"tz": "America\\/New_York", "gr": "MA", "g": "A6qOVH", "h": "wfLQtf", "l":
"orofrog", "al": "en-US,en;q=0.8", "hh": "1.usa.gov", "r":
"http:\\/\\/www.facebook.com\\/l\\/7AQEFzjSi\\/1.usa.gov\\/wfLQtf", "u":
"http:\\/\\/www.ncbi.nlm.nih.gov\\/pubmed\\/22415991", "t": 1331923247, "hc":
1331822918, "cy": "Danvers", "ll": [ 42.576698, -70.954903 ] }
Python 有内置和第三方库，用于将 JSON 字符串转换为 Python 字典。在这里，我们将使用json模块及其在我们下载的示例文件中的每一行上调用的loads函数：
import json
with open(path) as f:
records = [json.loads(line) for line in f]
结果对象records现在是一个 Python 字典列表：
In [18]: records[0]
Out[18]:
{'a': 'Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/535.11 (KHTML, like Gecko)
Chrome/17.0.963.78 Safari/535.11',
'al': 'en-US,en;q=0.8',
'c': 'US',
'cy': 'Danvers',
'g': 'A6qOVH',
'gr': 'MA',
'h': 'wfLQtf',
'hc': 1331822918,
'hh': '1.usa.gov',
'l': 'orofrog',
'll': [42.576698, -70.954903],
'nk': 1,
'r': 'http://www.facebook.com/l/7AQEFzjSi/1.usa.gov/wfLQtf',
't': 1331923247,
'tz': 'America/New_York',
'u': 'http://www.ncbi.nlm.nih.gov/pubmed/22415991'}
使用纯 Python 计算时区
假设我们有兴趣找出数据集中最常出现的时区（tz字段）。我们可以通过多种方式来实现这一点。首先，让我们再次使用列表推导式提取时区列表：
In [15]: time_zones = [rec["tz"] for rec in records]
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
<ipython-input-15-abdeba901c13> in <module>
----> 1 time_zones = [rec["tz"] for rec in records]
<ipython-input-15-abdeba901c13> in <listcomp>(.0)
----> 1 time_zones = [rec["tz"] for rec in records]
KeyError: 'tz'
糟糕！原来并非所有记录都有时区字段。我们可以通过在列表推导式末尾添加检查if "tz" in rec来处理这个问题：
In [16]: time_zones = [rec["tz"] for rec in records if "tz" in rec]
In [17]: time_zones[:10]
Out[17]:
['America/New_York',
'America/Denver',
'America/New_York',
'America/Sao_Paulo',
'America/New_York',
'America/New_York',
'Europe/Warsaw',
'',
'',
'']
仅查看前 10 个时区，我们会发现其中一些是未知的（空字符串）。您也可以将这些过滤掉，但我暂时保留它们。接下来，为了按时区生成计数，我将展示两种方法：一种更困难的方法（仅使用 Python 标准库）和一种更简单的方法（使用 pandas）。计数的一种方法是使用字典来存储计数，同时我们遍历时区：
def get_counts(sequence):
counts = {}
for x in sequence:
if x in counts:
counts[x] += 1
else:
counts[x] = 1
return counts
使用 Python 标准库中更高级的工具，您可以更简洁地编写相同的内容：
from collections import defaultdict
def get_counts2(sequence):
counts = defaultdict(int) # values will initialize to 0
for x in sequence:
counts[x] += 1
return counts
我将这个逻辑放在一个函数中，以使其更具可重用性。要在时区上使用它，只需传递time_zones列表：
In [20]: counts = get_counts(time_zones)
In [21]: counts["America/New_York"]
Out[21]: 1251
In [22]: len(time_zones)
Out[22]: 3440
如果我们想要前 10 个时区及其计数，我们可以通过(count, timezone)创建一个元组列表，并对其进行排序：
def top_counts(count_dict, n=10):
value_key_pairs = [(count, tz) for tz, count in count_dict.items()]
value_key_pairs.sort()
return value_key_pairs[-n:]
我们有：
In [24]: top_counts(counts)
Out[24]:
[(33, 'America/Sao_Paulo'),
(35, 'Europe/Madrid'),
(36, 'Pacific/Honolulu'),
(37, 'Asia/Tokyo'),
(74, 'Europe/London'),
(191, 'America/Denver'),
(382, 'America/Los_Angeles'),
(400, 'America/Chicago'),
(521, ''),
(1251, 'America/New_York')]
如果您搜索 Python 标准库，可能会找到collections.Counter类，这将使这个任务变得更简单：
In [25]: from collections import Counter
In [26]: counts = Counter(time_zones)
In [27]: counts.most_common(10)
Out[27]:
[('America/New_York', 1251),
('', 521),
('America/Chicago', 400),
('America/Los_Angeles', 382),
('America/Denver', 191),
('Europe/London', 74),
('Asia/Tokyo', 37),
('Pacific/Honolulu', 36),
('Europe/Madrid', 35),
('America/Sao_Paulo', 33)]
使用 pandas 计算时区
您可以通过将记录列表传递给pandas.DataFrame来从原始记录集创建一个 DataFrame：
In [28]: frame = pd.DataFrame(records)
我们可以查看有关这个新 DataFrame 的一些基本信息，比如列名、推断的列类型或缺失值的数量，使用frame.info()：
In [29]: frame.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 3560 entries, 0 to 3559
Data columns (total 18 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 a 3440 non-null object
1 c 2919 non-null object
2 nk 3440 non-null float64
3 tz 3440 non-null object
4 gr 2919 non-null object
5 g 3440 non-null object
6 h 3440 non-null object
7 l 3440 non-null object
8 al 3094 non-null object
9 hh 3440 non-null object
10 r 3440 non-null object
11 u 3440 non-null object
12 t 3440 non-null float64
13 hc 3440 non-null float64
14 cy 2919 non-null object
15 ll 2919 non-null object
16 _heartbeat_ 120 non-null float64
17 kw 93 non-null object
dtypes: float64(4), object(14)
memory usage: 500.8+ KB
In [30]: frame["tz"].head()
Out[30]:
0 America/New_York
1 America/Denver
2 America/New_York
3 America/Sao_Paulo
4 America/New_York
Name: tz, dtype: object
frame的输出显示为摘要视图，适用于大型 DataFrame 对象。然后我们可以使用 Series 的value_counts方法：
In [31]: tz_counts = frame["tz"].value_counts()
In [32]: tz_counts.head()
Out[32]:
tz
America/New_York 1251
521
America/Chicago 400
America/Los_Angeles 382
America/Denver 191
Name: count, dtype: int64
我们可以使用 matplotlib 可视化这些数据。我们可以通过为记录中的未知或缺失时区数据填充替代值来使图表更加美观。我们使用fillna方法替换缺失值，并使用布尔数组索引来处理空字符串：
In [33]: clean_tz = frame["tz"].fillna("Missing")
In [34]: clean_tz[clean_tz == ""] = "Unknown"
In [35]: tz_counts = clean_tz.value_counts()
In [36]: tz_counts.head()
Out[36]:
tz
America/New_York 1251
Unknown 521
America/Chicago 400
America/Los_Angeles 382
America/Denver 191
Name: count, dtype: int64
此时，我们可以使用seaborn 包制作一个水平条形图（参见 1.usa.gov 示例数据中的顶级时区以查看结果可视化）：
In [38]: import seaborn as sns
In [39]: subset = tz_counts.head()
In [40]: sns.barplot(y=subset.index, x=subset.to_numpy())
图 13.1：1.usa.gov 示例数据中的顶级时区
a字段包含有关用于执行 URL 缩短的浏览器、设备或应用程序的信息：
In [41]: frame["a"][1]
Out[41]: 'GoogleMaps/RochesterNY'
In [42]: frame["a"][50]
Out[42]: 'Mozilla/5.0 (Windows NT 5.1; rv:10.0.2) Gecko/20100101 Firefox/10.0.2'
In [43]: frame["a"][51][:50] # long line
Out[43]: 'Mozilla/5.0 (Linux; U; Android 2.2.2; en-us; LG-P9'
解析这些“代理”字符串中的所有有趣信息可能看起来是一项艰巨的任务。一种可能的策略是将字符串中的第一个标记（大致对应于浏览器功能）拆分出来，并对用户行为进行另一个摘要：
In [44]: results = pd.Series([x.split()[0] for x in frame["a"].dropna()])
In [45]: results.head(5)
Out[45]:
0 Mozilla/5.0
1 GoogleMaps/RochesterNY
2 Mozilla/4.0
3 Mozilla/5.0
4 Mozilla/5.0
dtype: object
In [46]: results.value_counts().head(8)
Out[46]:
Mozilla/5.0 2594
Mozilla/4.0 601
GoogleMaps/RochesterNY 121
Opera/9.80 34
TEST_INTERNET_AGENT 24
GoogleProducer 21
Mozilla/6.0 5
BlackBerry8520/5.0.0.681 4
Name: count, dtype: int64
现在，假设您想将顶级时区分解为 Windows 和非 Windows 用户。为简化起见，假设如果代理字符串中包含"Windows"字符串，则用户使用的是 Windows。由于一些代理缺失，我们将排除这些数据：
In [47]: cframe = frame[frame["a"].notna()].copy()
然后，我们想计算每行是否为 Windows 的值：
In [48]: cframe["os"] = np.where(cframe["a"].str.contains("Windows"),
....: "Windows", "Not Windows")
In [49]: cframe["os"].head(5)
Out[49]:
0 Windows
1 Not Windows
2 Windows
3 Not Windows
4 Windows
Name: os, dtype: object
然后，您可以按其时区列和这个新的操作系统列表对数据进行分组：
In [50]: by_tz_os = cframe.groupby(["tz", "os"])
类似于value_counts函数，可以使用size计算组计数。然后将结果重塑为表格，使用unstack：
In [51]: agg_counts = by_tz_os.size().unstack().fillna(0)
In [52]: agg_counts.head()
Out[52]:
os Not Windows Windows
tz
245.0 276.0
Africa/Cairo 0.0 3.0
Africa/Casablanca 0.0 1.0
Africa/Ceuta 0.0 2.0
Africa/Johannesburg 0.0 1.0
最后，让我们选择顶级的整体时区。为此，我从agg_counts中的行计数构建一个间接索引数组。在使用agg_counts.sum("columns")计算行计数后，我可以调用argsort()来获得一个可以用于升序排序的索引数组：
In [53]: indexer = agg_counts.sum("columns").argsort()
In [54]: indexer.values[:10]
Out[54]: array([24, 20, 21, 92, 87, 53, 54, 57, 26, 55])
我使用take按顺序选择行，然后切掉最后 10 行（最大值）：
In [55]: count_subset = agg_counts.take(indexer[-10:])
In [56]: count_subset
Out[56]:
os Not Windows Windows
tz
America/Sao_Paulo 13.0 20.0
Europe/Madrid 16.0 19.0
Pacific/Honolulu 0.0 36.0
Asia/Tokyo 2.0 35.0
Europe/London 43.0 31.0
America/Denver 132.0 59.0
America/Los_Angeles 130.0 252.0
America/Chicago 115.0 285.0
245.0 276.0
America/New_York 339.0 912.0
pandas 有一个方便的方法叫做nlargest，可以做同样的事情：
In [57]: agg_counts.sum(axis="columns").nlargest(10)
Out[57]:
tz
America/New_York 1251.0
521.0
America/Chicago 400.0
America/Los_Angeles 382.0
America/Denver 191.0
Europe/London 74.0
Asia/Tokyo 37.0
Pacific/Honolulu 36.0
Europe/Madrid 35.0
America/Sao_Paulo 33.0
dtype: float64
然后，可以绘制一个分组条形图，比较 Windows 和非 Windows 用户的数量，使用 seaborn 的barplot函数（参见按 Windows 和非 Windows 用户的顶级时区）。我首先调用count_subset.stack()并重置索引以重新排列数据，以便更好地与 seaborn 兼容：
In [59]: count_subset = count_subset.stack()
In [60]: count_subset.name = "total"
In [61]: count_subset = count_subset.reset_index()
In [62]: count_subset.head(10)
Out[62]:
tz os total
0 America/Sao_Paulo Not Windows 13.0
1 America/Sao_Paulo Windows 20.0
2 Europe/Madrid Not Windows 16.0
3 Europe/Madrid Windows 19.0
4 Pacific/Honolulu Not Windows 0.0
5 Pacific/Honolulu Windows 36.0
6 Asia/Tokyo Not Windows 2.0
7 Asia/Tokyo Windows 35.0
8 Europe/London Not Windows 43.0
9 Europe/London Windows 31.0
In [63]: sns.barplot(x="total", y="tz", hue="os", data=count_subset)
图 13.2：按 Windows 和非 Windows 用户的顶级时区
在较小的组中，很难看出 Windows 用户的相对百分比，因此让我们将组百分比归一化为 1：
def norm_total(group):
group["normed_total"] = group["total"] / group["total"].sum()
return group
results = count_subset.groupby("tz").apply(norm_total)
然后在出现频率最高的时区中 Windows 和非 Windows 用户的百分比中绘制这个图：
In [66]: sns.barplot(x="normed_total", y="tz", hue="os", data=results)
图 13.3：出现频率最高的时区中 Windows 和非 Windows 用户的百分比
我们可以通过使用transform方法和groupby更有效地计算归一化和：
In [67]: g = count_subset.groupby("tz")
In [68]: results2 = count_subset["total"] / g["total"].transform("sum")
13.2 MovieLens 1M 数据集
GroupLens Research提供了从 1990 年代末到 2000 年代初从 MovieLens 用户收集的多个电影评分数据集。数据提供了电影评分、电影元数据（类型和年份）以及关于用户的人口统计数据（年龄、邮政编码、性别认同和职业）。这些数据通常在基于机器学习算法的推荐系统的开发中很有兴趣。虽然我们在本书中没有详细探讨机器学习技术，但我将向您展示如何将这些数据集切分成您需要的确切形式。
MovieLens 1M 数据集包含从六千名用户对四千部电影收集的一百万个评分。它分布在三个表中：评分、用户信息和电影信息。我们可以使用pandas.read_table将每个表加载到一个 pandas DataFrame 对象中。在 Jupyter 单元格中运行以下代码：
unames = ["user_id", "gender", "age", "occupation", "zip"]
users = pd.read_table("datasets/movielens/users.dat", sep="::",
header=None, names=unames, engine="python")
rnames = ["user_id", "movie_id", "rating", "timestamp"]
ratings = pd.read_table("datasets/movielens/ratings.dat", sep="::",
header=None, names=rnames, engine="python")
mnames = ["movie_id", "title", "genres"]
movies = pd.read_table("datasets/movielens/movies.dat", sep="::",
header=None, names=mnames, engine="python")
您可以通过查看每个 DataFrame 来验证一切是否成功：
In [70]: users.head(5)
Out[70]:
user_id gender age occupation zip
0 1 F 1 10 48067
1 2 M 56 16 70072
2 3 M 25 15 55117
3 4 M 45 7 02460
4 5 M 25 20 55455
In [71]: ratings.head(5)
Out[71]:
user_id movie_id rating timestamp
0 1 1193 5 978300760
1 1 661 3 978302109
2 1 914 3 978301968
3 1 3408 4 978300275
4 1 2355 5 978824291
In [72]: movies.head(5)
Out[72]:
movie_id title genres
0 1 Toy Story (1995) Animation|Children's|Comedy
1 2 Jumanji (1995) Adventure|Children's|Fantasy
2 3 Grumpier Old Men (1995) Comedy|Romance
3 4 Waiting to Exhale (1995) Comedy|Drama
4 5 Father of the Bride Part II (1995) Comedy
In [73]: ratings
Out[73]:
user_id movie_id rating timestamp
0 1 1193 5 978300760
1 1 661 3 978302109
2 1 914 3 978301968
3 1 3408 4 978300275
4 1 2355 5 978824291
... ... ... ... ...
1000204 6040 1091 1 956716541
1000205 6040 1094 5 956704887
1000206 6040 562 5 956704746
1000207 6040 1096 4 956715648
1000208 6040 1097 4 956715569
[1000209 rows x 4 columns]
请注意，年龄和职业被编码为整数，表示数据集的README文件中描述的组。分析分布在三个表中的数据并不是一项简单的任务；例如，假设您想要按性别身份和年龄计算特定电影的平均评分。正如您将看到的，将所有数据合并到一个单一表中更方便。使用 pandas 的merge函数，我们首先将ratings与users合并，然后将该结果与movies数据合并。pandas 根据重叠的名称推断要用作合并（或join）键的列：
In [74]: data = pd.merge(pd.merge(ratings, users), movies)
In [75]: data
Out[75]:
user_id movie_id rating timestamp gender age occupation zip
0 1 1193 5 978300760 F 1 10 48067 \
1 2 1193 5 978298413 M 56 16 70072
2 12 1193 4 978220179 M 25 12 32793
3 15 1193 4 978199279 M 25 7 22903
4 17 1193 5 978158471 M 50 1 95350
... ... ... ... ... ... ... ... ...
1000204 5949 2198 5 958846401 M 18 17 47901
1000205 5675 2703 3 976029116 M 35 14 30030
1000206 5780 2845 1 958153068 M 18 17 92886
1000207 5851 3607 5 957756608 F 18 20 55410
1000208 5938 2909 4 957273353 M 25 1 35401
title genres
0 One Flew Over the Cuckoo's Nest (1975) Drama
1 One Flew Over the Cuckoo's Nest (1975) Drama
2 One Flew Over the Cuckoo's Nest (1975) Drama
3 One Flew Over the Cuckoo's Nest (1975) Drama
4 One Flew Over the Cuckoo's Nest (1975) Drama
... ... ...
1000204 Modulations (1998) Documentary
1000205 Broken Vessels (1998) Drama
1000206 White Boys (1999) Drama
1000207 One Little Indian (1973) Comedy|Drama|Western
1000208 Five Wives, Three Secretaries and Me (1998) Documentary
[1000209 rows x 10 columns]
In [76]: data.iloc[0]
Out[76]:
user_id 1
movie_id 1193
rating 5
timestamp 978300760
gender F
age 1
occupation 10
zip 48067
title One Flew Over the Cuckoo's Nest (1975)
genres Drama
Name: 0, dtype: object
为了获得按性别分组的每部电影的平均评分，我们可以使用pivot_table方法：
In [77]: mean_ratings = data.pivot_table("rating", index="title",
....: columns="gender", aggfunc="mean")
In [78]: mean_ratings.head(5)
Out[78]:
gender F M
title
$1,000,000 Duck (1971) 3.375000 2.761905
'Night Mother (1986) 3.388889 3.352941
'Til There Was You (1997) 2.675676 2.733333
'burbs, The (1989) 2.793478 2.962085
...And Justice for All (1979) 3.828571 3.689024
这产生了另一个包含平均评分的 DataFrame，其中电影标题作为行标签（“索引”），性别作为列标签。我首先筛选出至少收到 250 个评分的电影（一个任意的数字）；为此，我按标题对数据进行分组，并使用size()来获取每个标题的组大小的 Series：
In [79]: ratings_by_title = data.groupby("title").size()
In [80]: ratings_by_title.head()
Out[80]:
title
$1,000,000 Duck (1971) 37
'Night Mother (1986) 70
'Til There Was You (1997) 52
'burbs, The (1989) 303
...And Justice for All (1979) 199
dtype: int64
In [81]: active_titles = ratings_by_title.index[ratings_by_title >= 250]
In [82]: active_titles
Out[82]:
Index([''burbs, The (1989)', '10 Things I Hate About You (1999)',
'101 Dalmatians (1961)', '101 Dalmatians (1996)', '12 Angry Men (1957)',
'13th Warrior, The (1999)', '2 Days in the Valley (1996)',
'20,000 Leagues Under the Sea (1954)', '2001: A Space Odyssey (1968)',
'2010 (1984)',
...
'X-Men (2000)', 'Year of Living Dangerously (1982)',
'Yellow Submarine (1968)', 'You've Got Mail (1998)',
'Young Frankenstein (1974)', 'Young Guns (1988)',
'Young Guns II (1990)', 'Young Sherlock Holmes (1985)',
'Zero Effect (1998)', 'eXistenZ (1999)'],
dtype='object', name='title', length=1216)
然后，可以使用至少收到 250 个评分的标题的索引来从mean_ratings中选择行，使用.loc：
In [83]: mean_ratings = mean_ratings.loc[active_titles]
In [84]: mean_ratings
Out[84]:
gender F M
title
'burbs, The (1989) 2.793478 2.962085
10 Things I Hate About You (1999) 3.646552 3.311966
101 Dalmatians (1961) 3.791444 3.500000
101 Dalmatians (1996) 3.240000 2.911215
12 Angry Men (1957) 4.184397 4.328421
... ... ...
Young Guns (1988) 3.371795 3.425620
Young Guns II (1990) 2.934783 2.904025
Young Sherlock Holmes (1985) 3.514706 3.363344
Zero Effect (1998) 3.864407 3.723140
eXistenZ (1999) 3.098592 3.289086
[1216 rows x 2 columns]
要查看女性观众最喜欢的电影，我们可以按降序排序F列：
In [86]: top_female_ratings = mean_ratings.sort_values("F", ascending=False)
In [87]: top_female_ratings.head()
Out[87]:
gender F M
title
Close Shave, A (1995) 4.644444 4.473795
Wrong Trousers, The (1993) 4.588235 4.478261
Sunset Blvd. (a.k.a. Sunset Boulevard) (1950) 4.572650 4.464589
Wallace & Gromit: The Best of Aardman Animation (1996) 4.563107 4.385075
Schindler's List (1993) 4.562602 4.491415
测量评分分歧
假设您想要找到在男性和女性观众之间最具分歧的电影。一种方法是向mean_ratings添加一个包含平均值差异的列，然后按照该列进行排序：
In [88]: mean_ratings["diff"] = mean_ratings["M"] - mean_ratings["F"]
按照"diff"排序，可以得到评分差异最大的电影，以便看到哪些电影更受女性喜欢：
In [89]: sorted_by_diff = mean_ratings.sort_values("diff")
In [90]: sorted_by_diff.head()
Out[90]:
gender F M diff
title
Dirty Dancing (1987) 3.790378 2.959596 -0.830782
Jumpin' Jack Flash (1986) 3.254717 2.578358 -0.676359
Grease (1978) 3.975265 3.367041 -0.608224
Little Women (1994) 3.870588 3.321739 -0.548849
Steel Magnolias (1989) 3.901734 3.365957 -0.535777
颠倒行的顺序并再次切片前 10 行，我们得到了男性喜欢但女性评分不高的电影：
In [91]: sorted_by_diff[::-1].head()
Out[91]:
gender F M diff
title
Good, The Bad and The Ugly, The (1966) 3.494949 4.221300 0.726351
Kentucky Fried Movie, The (1977) 2.878788 3.555147 0.676359
Dumb & Dumber (1994) 2.697987 3.336595 0.638608
Longest Day, The (1962) 3.411765 4.031447 0.619682
Cable Guy, The (1996) 2.250000 2.863787 0.613787
假设您想要找到在观众中引起最大分歧的电影，而不考虑性别认同。分歧可以通过评分的方差或标准差来衡量。为了得到这个结果，我们首先按标题计算评分的标准差，然后筛选出活跃的标题：
In [92]: rating_std_by_title = data.groupby("title")["rating"].std()
In [93]: rating_std_by_title = rating_std_by_title.loc[active_titles]
In [94]: rating_std_by_title.head()
Out[94]:
title
'burbs, The (1989) 1.107760
10 Things I Hate About You (1999) 0.989815
101 Dalmatians (1961) 0.982103
101 Dalmatians (1996) 1.098717
12 Angry Men (1957) 0.812731
Name: rating, dtype: float64
然后，我们按降序排序并选择前 10 行，这大致是评分最具分歧的 10 部电影：
In [95]: rating_std_by_title.sort_values(ascending=False)[:10]
Out[95]:
title
Dumb & Dumber (1994) 1.321333
Blair Witch Project, The (1999) 1.316368
Natural Born Killers (1994) 1.307198
Tank Girl (1995) 1.277695
Rocky Horror Picture Show, The (1975) 1.260177
Eyes Wide Shut (1999) 1.259624
Evita (1996) 1.253631
Billy Madison (1995) 1.249970
Fear and Loathing in Las Vegas (1998) 1.246408
Bicentennial Man (1999) 1.245533
Name: rating, dtype: float64
您可能已经注意到电影类型是以管道分隔（|）的字符串给出的，因为一部电影可以属于多种类型。为了帮助我们按类型对评分数据进行分组，我们可以在 DataFrame 上使用explode方法。让我们看看这是如何工作的。首先，我们可以使用 Series 上的str.split方法将类型字符串拆分为类型列表：
In [96]: movies["genres"].head()
Out[96]:
0 Animation|Children's|Comedy
1 Adventure|Children's|Fantasy
2 Comedy|Romance
3 Comedy|Drama
4 Comedy
Name: genres, dtype: object
In [97]: movies["genres"].head().str.split("|")
Out[97]:
0 [Animation, Children's, Comedy]
1 [Adventure, Children's, Fantasy]
2 [Comedy, Romance]
3 [Comedy, Drama]
4 [Comedy]
Name: genres, dtype: object
In [98]: movies["genre"] = movies.pop("genres").str.split("|")
In [99]: movies.head()
Out[99]:
movie_id title
0 1 Toy Story (1995) \
1 2 Jumanji (1995)
2 3 Grumpier Old Men (1995)
3 4 Waiting to Exhale (1995)
4 5 Father of the Bride Part II (1995)
genre
0 [Animation, Children's, Comedy]
1 [Adventure, Children's, Fantasy]
2 [Comedy, Romance]
3 [Comedy, Drama]
4 [Comedy]
现在，调用movies.explode("genre")会生成一个新的 DataFrame，其中每个电影类型列表中的“内部”元素都有一行。例如，如果一部电影被分类为喜剧和浪漫片，那么结果中将有两行，一行只有“喜剧”，另一行只有“浪漫片”：
In [100]: movies_exploded = movies.explode("genre")
In [101]: movies_exploded[:10]
Out[101]:
movie_id title genre
0 1 Toy Story (1995) Animation
0 1 Toy Story (1995) Children's
0 1 Toy Story (1995) Comedy
1 2 Jumanji (1995) Adventure
1 2 Jumanji (1995) Children's
1 2 Jumanji (1995) Fantasy
2 3 Grumpier Old Men (1995) Comedy
2 3 Grumpier Old Men (1995) Romance
3 4 Waiting to Exhale (1995) Comedy
3 4 Waiting to Exhale (1995) Drama
现在，我们可以将所有三个表合并在一起，并按类型分组：
In [102]: ratings_with_genre = pd.merge(pd.merge(movies_exploded, ratings), users
)
In [103]: ratings_with_genre.iloc[0]
Out[103]:
movie_id 1
title Toy Story (1995)
genre Animation
user_id 1
rating 5
timestamp 978824268
gender F
age 1
occupation 10
zip 48067
Name: 0, dtype: object
In [104]: genre_ratings = (ratings_with_genre.groupby(["genre", "age"])
.....: ["rating"].mean()
.....: .unstack("age"))
In [105]: genre_ratings[:10]
Out[105]:
age 1 18 25 35 45 50
genre
Action 3.506385 3.447097 3.453358 3.538107 3.528543 3.611333 \
Adventure 3.449975 3.408525 3.443163 3.515291 3.528963 3.628163
Animation 3.476113 3.624014 3.701228 3.740545 3.734856 3.780020
Children's 3.241642 3.294257 3.426873 3.518423 3.527593 3.556555
Comedy 3.497491 3.460417 3.490385 3.561984 3.591789 3.646868
Crime 3.710170 3.668054 3.680321 3.733736 3.750661 3.810688
Documentary 3.730769 3.865865 3.946690 3.953747 3.966521 3.908108
Drama 3.794735 3.721930 3.726428 3.782512 3.784356 3.878415
Fantasy 3.317647 3.353778 3.452484 3.482301 3.532468 3.581570
Film-Noir 4.145455 3.997368 4.058725 4.064910 4.105376 4.175401
age 56
genre
Action 3.610709
Adventure 3.649064
Animation 3.756233
Children's 3.621822
Comedy 3.650949
Crime 3.832549
Documentary 3.961538
Drama 3.933465
Fantasy 3.532700
Film-Noir 4.125932
13.3 美国婴儿姓名 1880-2010
美国社会保障管理局（SSA）提供了从 1880 年到现在的婴儿名字频率数据。Hadley Wickham，几个流行 R 包的作者，在 R 中说明数据操作时使用了这个数据集。
我们需要进行一些数据整理来加载这个数据集，但一旦我们这样做了，我们将得到一个看起来像这样的 DataFrame：
In [4]: names.head(10)
Out[4]:
name sex births year
0 Mary F 7065 1880
1 Anna F 2604 1880
2 Emma F 2003 1880
3 Elizabeth F 1939 1880
4 Minnie F 1746 1880
5 Margaret F 1578 1880
6 Ida F 1472 1880
7 Alice F 1414 1880
8 Bertha F 1320 1880
9 Sarah F 1288 1880
有许多事情你可能想要对数据集做：
-
可视化随着时间推移给定名字（您自己的名字或其他名字）的婴儿比例
-
确定一个名字的相对排名
-
确定每年最受欢迎的名字或受欢迎程度增长或下降最多的名字
-
分析名字的趋势：元音、辅音、长度、整体多样性、拼写变化、首尾字母
-
分析趋势的外部来源：圣经名字、名人、人口统计学
使用本书中的工具，许多这类分析都可以实现，所以我会带你走一些。
截至目前，美国社会保障管理局提供了数据文件，每年一个文件，其中包含每个性别/名字组合的总出生数。您可以下载这些文件的原始存档。
如果您在阅读此页面时发现已移动，很可能可以通过互联网搜索再次找到。下载“国家数据”文件names.zip并解压缩后，您将获得一个包含一系列文件如yob1880.txt的目录。我使用 Unix 的head命令查看其中一个文件的前 10 行（在 Windows 上，您可以使用more命令或在文本编辑器中打开）：
In [106]: !head -n 10 datasets/babynames/yob1880.txt
Mary,F,7065
Anna,F,2604
Emma,F,2003
Elizabeth,F,1939
Minnie,F,1746
Margaret,F,1578
Ida,F,1472
Alice,F,1414
Bertha,F,1320
Sarah,F,1288
由于这已经是逗号分隔形式，可以使用pandas.read_csv将其加载到 DataFrame 中：
In [107]: names1880 = pd.read_csv("datasets/babynames/yob1880.txt",
.....: names=["name", "sex", "births"])
In [108]: names1880
Out[108]:
name sex births
0 Mary F 7065
1 Anna F 2604
2 Emma F 2003
3 Elizabeth F 1939
4 Minnie F 1746
... ... .. ...
1995 Woodie M 5
1996 Worthy M 5
1997 Wright M 5
1998 York M 5
1999 Zachariah M 5
[2000 rows x 3 columns]
这些文件只包含每年至少有五次出现的名字，所以为了简单起见，我们可以使用按性别的出生列的总和作为该年出生的总数：
In [109]: names1880.groupby("sex")["births"].sum()
Out[109]:
sex
F 90993
M 110493
Name: births, dtype: int64
由于数据集按年份分成文件，首先要做的事情之一是将所有数据组装到一个单独的 DataFrame 中，并进一步添加一个year字段。您可以使用pandas.concat来做到这一点。在 Jupyter 单元格中运行以下内容：
pieces = []
for year in range(1880, 2011):
path = f"datasets/babynames/yob{year}.txt"
frame = pd.read_csv(path, names=["name", "sex", "births"])
# Add a column for the year
frame["year"] = year
pieces.append(frame)
# Concatenate everything into a single DataFrame
names = pd.concat(pieces, ignore_index=True)
这里有几件事情需要注意。首先，记住concat默认按行组合 DataFrame 对象。其次，您必须传递ignore_index=True，因为我们不关心从pandas.read_csv返回的原始行号。因此，现在我们有一个包含所有年份的所有名字数据的单个 DataFrame：
In [111]: names
Out[111]:
name sex births year
0 Mary F 7065 1880
1 Anna F 2604 1880
2 Emma F 2003 1880
3 Elizabeth F 1939 1880
4 Minnie F 1746 1880
... ... .. ... ...
1690779 Zymaire M 5 2010
1690780 Zyonne M 5 2010
1690781 Zyquarius M 5 2010
1690782 Zyran M 5 2010
1690783 Zzyzx M 5 2010
[1690784 rows x 4 columns]
有了这些数据，我们可以开始使用groupby或pivot_table在年份和性别水平上对数据进行聚合（参见按性别和年份统计的总出生数）：
In [112]: total_births = names.pivot_table("births", index="year",
.....: columns="sex", aggfunc=sum)
In [113]: total_births.tail()
Out[113]:
sex F M
year
2006 1896468 2050234
2007 1916888 2069242
2008 1883645 2032310
2009 1827643 1973359
2010 1759010 1898382
In [114]: total_births.plot(title="Total births by sex and year")
图 13.4：按性别和年份统计的总出生数
接下来，让我们插入一个名为prop的列，该列显示每个名字相对于总出生数的比例。prop值为0.02表示每 100 个婴儿中有 2 个被赋予特定的名字。因此，我们按年份和性别对数据进行分组，然后向每个组添加新列：
def add_prop(group):
group["prop"] = group["births"] / group["births"].sum()
return group
names = names.groupby(["year", "sex"], group_keys=False).apply(add_prop)
现在得到的完整数据集现在具有以下列：
In [116]: names
Out[116]:
name sex births year prop
0 Mary F 7065 1880 0.077643
1 Anna F 2604 1880 0.028618
2 Emma F 2003 1880 0.022013
3 Elizabeth F 1939 1880 0.021309
4 Minnie F 1746 1880 0.019188
... ... .. ... ... ...
1690779 Zymaire M 5 2010 0.000003
1690780 Zyonne M 5 2010 0.000003
1690781 Zyquarius M 5 2010 0.000003
1690782 Zyran M 5 2010 0.000003
1690783 Zzyzx M 5 2010 0.000003
[1690784 rows x 5 columns]
在执行这样的组操作时，通常很有价值进行一些合理性检查，比如验证所有组中prop列的总和是否为 1：
In [117]: names.groupby(["year", "sex"])["prop"].sum()
Out[117]:
year sex
1880 F 1.0
M 1.0
1881 F 1.0
M 1.0
1882 F 1.0
...
2008 M 1.0
2009 F 1.0
M 1.0
2010 F 1.0
M 1.0
Name: prop, Length: 262, dtype: float64
现在这样做了，我将提取数据的一个子集以便进一步分析：每个性别/年份组合的前 1000 个名字。这是另一个组操作：
In [118]: def get_top1000(group):
.....: return group.sort_values("births", ascending=False)[:1000]
In [119]: grouped = names.groupby(["year", "sex"])
In [120]: top1000 = grouped.apply(get_top1000)
In [121]: top1000.head()
Out[121]:
name sex births year prop
year sex
1880 F 0 Mary F 7065 1880 0.077643
1 Anna F 2604 1880 0.028618
2 Emma F 2003 1880 0.022013
3 Elizabeth F 1939 1880 0.021309
4 Minnie F 1746 1880 0.019188
我们可以删除组索引，因为我们不需要它进行分析：
In [122]: top1000 = top1000.reset_index(drop=True)
现在得到的数据集要小得多：
In [123]: top1000.head()
Out[123]:
name sex births year prop
0 Mary F 7065 1880 0.077643
1 Anna F 2604 1880 0.028618
2 Emma F 2003 1880 0.022013
3 Elizabeth F 1939 1880 0.021309
4 Minnie F 1746 1880 0.019188
我们将在接下来的数据调查中使用这个前一千个数据集。
分析命名趋势
有了完整的数据集和前一千个数据集，我们可以开始分析各种有趣的命名趋势。首先，我们可以将前一千个名字分为男孩和女孩部分：
In [124]: boys = top1000[top1000["sex"] == "M"]
In [125]: girls = top1000[top1000["sex"] == "F"]
简单的时间序列，比如每年约翰或玛丽的数量，可以绘制，但需要一些操作才能更有用。让我们形成一个按年份和姓名总数的数据透视表：
In [126]: total_births = top1000.pivot_table("births", index="year",
.....: columns="name",
.....: aggfunc=sum)
现在，可以使用 DataFrame 的plot方法为一些名字绘制图表（一些男孩和女孩名字随时间变化显示了结果）：
In [127]: total_births.info()
<class 'pandas.core.frame.DataFrame'>
Index: 131 entries, 1880 to 2010
Columns: 6868 entries, Aaden to Zuri
dtypes: float64(6868)
memory usage: 6.9 MB
In [128]: subset = total_births[["John", "Harry", "Mary", "Marilyn"]]
In [129]: subset.plot(subplots=True, figsize=(12, 10),
.....: title="Number of births per year")
图 13.5：一些男孩和女孩名字随时间变化
看到这个，你可能会得出结论，这些名字已经不再受到美国人口的青睐。但事实实际上比这更复杂，将在下一节中探讨。
衡量命名多样性的增加
减少图表的原因之一是越来越少的父母选择常见的名字给他们的孩子。这个假设可以在数据中进行探索和确认。一个度量是由前 1000 个最受欢迎的名字代表的出生比例，我按年份和性别进行汇总和绘制（性别在前一千个名字中所代表的出生比例显示了结果图）：
In [131]: table = top1000.pivot_table("prop", index="year",
.....: columns="sex", aggfunc=sum)
In [132]: table.plot(title="Sum of table1000.prop by year and sex",
.....: yticks=np.linspace(0, 1.2, 13))
图 13.6：性别在前一千个名字中所代表的出生比例
您可以看到，确实存在着越来越多的名字多样性（前一千名中总比例减少）。另一个有趣的指标是在出生的前 50%中按照从高到低的流行度顺序取的不同名字的数量。这个数字更难计算。让我们只考虑 2010 年的男孩名字：
In [133]: df = boys[boys["year"] == 2010]
In [134]: df
Out[134]:
name sex births year prop
260877 Jacob M 21875 2010 0.011523
260878 Ethan M 17866 2010 0.009411
260879 Michael M 17133 2010 0.009025
260880 Jayden M 17030 2010 0.008971
260881 William M 16870 2010 0.008887
... ... .. ... ... ...
261872 Camilo M 194 2010 0.000102
261873 Destin M 194 2010 0.000102
261874 Jaquan M 194 2010 0.000102
261875 Jaydan M 194 2010 0.000102
261876 Maxton M 193 2010 0.000102
[1000 rows x 5 columns]
在对prop进行降序排序后，我们想知道最受欢迎的名字中有多少个名字达到了 50%。您可以编写一个for循环来执行此操作，但使用矢量化的 NumPy 方法更具计算效率。对prop进行累积求和cumsum，然后调用searchsorted方法返回0.5需要插入的累积和位置，以保持其按顺序排序：
In [135]: prop_cumsum = df["prop"].sort_values(ascending=False).cumsum()
In [136]: prop_cumsum[:10]
Out[136]:
260877 0.011523
260878 0.020934
260879 0.029959
260880 0.038930
260881 0.047817
260882 0.056579
260883 0.065155
260884 0.073414
260885 0.081528
260886 0.089621
Name: prop, dtype: float64
In [137]: prop_cumsum.searchsorted(0.5)
Out[137]: 116
由于数组是从零开始索引的，将此结果加 1 将得到 117 的结果。相比之下，在 1900 年，这个数字要小得多：
In [138]: df = boys[boys.year == 1900]
In [139]: in1900 = df.sort_values("prop", ascending=False).prop.cumsum()
In [140]: in1900.searchsorted(0.5) + 1
Out[140]: 25
现在，您可以将此操作应用于每个年份/性别组合，对这些字段进行groupby，并apply一个返回每个组计数的函数：
def get_quantile_count(group, q=0.5):
group = group.sort_values("prop", ascending=False)
return group.prop.cumsum().searchsorted(q) + 1
diversity = top1000.groupby(["year", "sex"]).apply(get_quantile_count)
diversity = diversity.unstack()
这个结果 DataFrame diversity 现在有两个时间序列，一个用于每个性别，按年份索引。这可以像以前一样进行检查和绘制（参见按年份绘制的多样性指标）：
In [143]: diversity.head()
Out[143]:
sex F M
year
1880 38 14
1881 38 14
1882 38 15
1883 39 15
1884 39 16
In [144]: diversity.plot(title="Number of popular names in top 50%")
图 13.7：按年份绘制的多样性指标
正如你所看到的，女孩名字一直比男孩名字更多样化，而且随着时间的推移，它们变得更加多样化。关于到底是什么推动了这种多样性的进一步分析，比如替代拼写的增加，留给读者自行探讨。
“最后一个字母”革命
在 2007 年，婴儿姓名研究员劳拉·瓦滕伯格指出，过去 100 年来，以最后一个字母结尾的男孩名字的分布发生了显著变化。为了看到这一点，我们首先按年份、性别和最后一个字母聚合完整数据集中的所有出生情况：
def get_last_letter(x):
return x[-1]
last_letters = names["name"].map(get_last_letter)
last_letters.name = "last_letter"
table = names.pivot_table("births", index=last_letters,
columns=["sex", "year"], aggfunc=sum)
然后我们选择三个代表性年份跨越历史，并打印前几行：
In [146]: subtable = table.reindex(columns=[1910, 1960, 2010], level="year")
In [147]: subtable.head()
Out[147]:
sex F M
year 1910 1960 2010 1910 1960 2010
last_letter
a 108376.0 691247.0 670605.0 977.0 5204.0 28438.0
b NaN 694.0 450.0 411.0 3912.0 38859.0
c 5.0 49.0 946.0 482.0 15476.0 23125.0
d 6750.0 3729.0 2607.0 22111.0 262112.0 44398.0
e 133569.0 435013.0 313833.0 28655.0 178823.0 129012.0
接下来，通过总出生数对表进行标准化，计算一个包含每个性别以每个字母结尾的总出生比例的新表：
In [148]: subtable.sum()
Out[148]:
sex year
F 1910 396416.0
1960 2022062.0
2010 1759010.0
M 1910 194198.0
1960 2132588.0
2010 1898382.0
dtype: float64
In [149]: letter_prop = subtable / subtable.sum()
In [150]: letter_prop
Out[150]:
sex F M
year 1910 1960 2010 1910 1960 2010
last_letter
a 0.273390 0.341853 0.381240 0.005031 0.002440 0.014980
b NaN 0.000343 0.000256 0.002116 0.001834 0.020470
c 0.000013 0.000024 0.000538 0.002482 0.007257 0.012181
d 0.017028 0.001844 0.001482 0.113858 0.122908 0.023387
e 0.336941 0.215133 0.178415 0.147556 0.083853 0.067959
... ... ... ... ... ... ...
v NaN 0.000060 0.000117 0.000113 0.000037 0.001434
w 0.000020 0.000031 0.001182 0.006329 0.007711 0.016148
x 0.000015 0.000037 0.000727 0.003965 0.001851 0.008614
y 0.110972 0.152569 0.116828 0.077349 0.160987 0.058168
z 0.002439 0.000659 0.000704 0.000170 0.000184 0.001831
[26 rows x 6 columns]
现在有了字母比例，我们可以按年份将每个性别分解为条形图（参见以每个字母结尾的男孩和女孩名字的比例）：
import matplotlib.pyplot as plt
fig, axes = plt.subplots(2, 1, figsize=(10, 8))
letter_prop["M"].plot(kind="bar", rot=0, ax=axes[0], title="Male")
letter_prop["F"].plot(kind="bar", rot=0, ax=axes[1], title="Female",
legend=False)
图 13.8：以每个字母结尾的男孩和女孩名字的比例
正如您所看到的，自 20 世纪 60 年代以来，以n结尾的男孩名字经历了显著增长。回到之前创建的完整表格，再次按年份和性别进行标准化，并选择男孩名字的一部分字母，最后转置使每一列成为一个时间序列：
In [153]: letter_prop = table / table.sum()
In [154]: dny_ts = letter_prop.loc[["d", "n", "y"], "M"].T
In [155]: dny_ts.head()
Out[155]:
last_letter d n y
year
1880 0.083055 0.153213 0.075760
1881 0.083247 0.153214 0.077451
1882 0.085340 0.149560 0.077537
1883 0.084066 0.151646 0.079144
1884 0.086120 0.149915 0.080405
有了这个时间序列的 DataFrame，我可以再次使用其plot方法制作时间趋势图（请参见随时间变化以 d/n/y 结尾的男孩出生比例）：
In [158]: dny_ts.plot()
图 13.9：随时间变化以 d/n/y 结尾的男孩出生比例
男孩名字变成女孩名字（反之亦然）
另一个有趣的趋势是查看在样本早期更受一性别欢迎，但随着时间推移已成为另一性别的首选名字的名字。一个例子是 Lesley 或 Leslie 这个名字。回到top1000 DataFrame，我计算出数据集中以"Lesl"开头的名字列表：
In [159]: all_names = pd.Series(top1000["name"].unique())
In [160]: lesley_like = all_names[all_names.str.contains("Lesl")]
In [161]: lesley_like
Out[161]:
632 Leslie
2294 Lesley
4262 Leslee
4728 Lesli
6103 Lesly
dtype: object
然后，我们可以筛选出那些名字，按名字分组对出生进行求和，以查看相对频率：
In [162]: filtered = top1000[top1000["name"].isin(lesley_like)]
In [163]: filtered.groupby("name")["births"].sum()
Out[163]:
name
Leslee 1082
Lesley 35022
Lesli 929
Leslie 370429
Lesly 10067
Name: births, dtype: int64
接下来，让我们按性别和年份进行聚合，并在年份内进行归一化：
In [164]: table = filtered.pivot_table("births", index="year",
.....: columns="sex", aggfunc="sum")
In [165]: table = table.div(table.sum(axis="columns"), axis="index")
In [166]: table.tail()
Out[166]:
sex F M
year
2006 1.0 NaN
2007 1.0 NaN
2008 1.0 NaN
2009 1.0 NaN
2010 1.0 NaN
最后，现在可以制作按性别随时间变化的分布图（请参见随时间变化男/女 Lesley 样式名字的比例）：
In [168]: table.plot(style={"M": "k-", "F": "k--"})
图 13.10：随时间变化男/女 Lesley 样式名字的比例
13.4 USDA 食品数据库
美国农业部（USDA）提供了一个食品营养信息数据库。程序员 Ashley Williams 以 JSON 格式创建了这个数据库的一个版本。记录看起来像这样：
{
"id": 21441,
"description": "KENTUCKY FRIED CHICKEN, Fried Chicken, EXTRA CRISPY,
Wing, meat and skin with breading",
"tags": ["KFC"],
"manufacturer": "Kentucky Fried Chicken",
"group": "Fast Foods",
"portions": [
{
"amount": 1,
"unit": "wing, with skin",
"grams": 68.0
},
...
],
"nutrients": [
{
"value": 20.8,
"units": "g",
"description": "Protein",
"group": "Composition"
},
...
]
}
每种食物都有一些标识属性，还有两个营养素和分量大小的列表。这种形式的数据不太适合分析，因此我们需要做一些工作，将数据整理成更好的形式。
您可以使用您选择的任何 JSON 库将此文件加载到 Python 中。我将使用内置的 Python json模块：
In [169]: import json
In [170]: db = json.load(open("datasets/usda_food/database.json"))
In [171]: len(db)
Out[171]: 6636
db中的每个条目都是一个包含单个食物所有数据的字典。"nutrients"字段是一个字典列表，每个营养素一个：
In [172]: db[0].keys()
Out[172]: dict_keys(['id', 'description', 'tags', 'manufacturer', 'group', 'porti
ons', 'nutrients'])
In [173]: db[0]["nutrients"][0]
Out[173]:
{'value': 25.18,
'units': 'g',
'description': 'Protein',
'group': 'Composition'}
In [174]: nutrients = pd.DataFrame(db[0]["nutrients"])
In [175]: nutrients.head(7)
Out[175]:
value units description group
0 25.18 g Protein Composition
1 29.20 g Total lipid (fat) Composition
2 3.06 g Carbohydrate, by difference Composition
3 3.28 g Ash Other
4 376.00 kcal Energy Energy
5 39.28 g Water Composition
6 1573.00 kJ Energy Energy
将字典列表转换为 DataFrame 时，我们可以指定要提取的字段列表。我们将提取食物名称、组、ID 和制造商：
In [176]: info_keys = ["description", "group", "id", "manufacturer"]
In [177]: info = pd.DataFrame(db, columns=info_keys)
In [178]: info.head()
Out[178]:
description group id
0 Cheese, caraway Dairy and Egg Products 1008 \
1 Cheese, cheddar Dairy and Egg Products 1009
2 Cheese, edam Dairy and Egg Products 1018
3 Cheese, feta Dairy and Egg Products 1019
4 Cheese, mozzarella, part skim milk Dairy and Egg Products 1028
manufacturer
0
1
2
3
4
In [179]: info.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 6636 entries, 0 to 6635
Data columns (total 4 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 description 6636 non-null object
1 group 6636 non-null object
2 id 6636 non-null int64
3 manufacturer 5195 non-null object
dtypes: int64(1), object(3)
memory usage: 207.5+ KB
从info.info()的输出中，我们可以看到manufacturer列中有缺失数据。
您可以使用value_counts查看食物组的分布：
In [180]: pd.value_counts(info["group"])[:10]
Out[180]:
group
Vegetables and Vegetable Products 812
Beef Products 618
Baked Products 496
Breakfast Cereals 403
Legumes and Legume Products 365
Fast Foods 365
Lamb, Veal, and Game Products 345
Sweets 341
Fruits and Fruit Juices 328
Pork Products 328
Name: count, dtype: int64
现在，要对所有营养数据进行一些分析，最简单的方法是将每种食物的营养成分组装成一个单独的大表格。为此，我们需要采取几个步骤。首先，我将把每个食物营养列表转换为一个 DataFrame，添加一个食物id的列，并将 DataFrame 附加到列表中。然后，可以使用concat将它们连接起来。在 Jupyter 单元格中运行以下代码：
nutrients = []
for rec in db:
fnuts = pd.DataFrame(rec["nutrients"])
fnuts["id"] = rec["id"]
nutrients.append(fnuts)
nutrients = pd.concat(nutrients, ignore_index=True)
如果一切顺利，nutrients应该是这样的：
In [182]: nutrients
Out[182]:
value units description group id
0 25.180 g Protein Composition 1008
1 29.200 g Total lipid (fat) Composition 1008
2 3.060 g Carbohydrate, by difference Composition 1008
3 3.280 g Ash Other 1008
4 376.000 kcal Energy Energy 1008
... ... ... ... ... ...
389350 0.000 mcg Vitamin B-12, added Vitamins 43546
389351 0.000 mg Cholesterol Other 43546
389352 0.072 g Fatty acids, total saturated Other 43546
389353 0.028 g Fatty acids, total monounsaturated Other 43546
389354 0.041 g Fatty acids, total polyunsaturated Other 43546
[389355 rows x 5 columns]
我注意到这个 DataFrame 中有重复项，所以删除它们会更容易：
In [183]: nutrients.duplicated().sum() # number of duplicates
Out[183]: 14179
In [184]: nutrients = nutrients.drop_duplicates()
由于 DataFrame 对象中都有"group"和"description"，我们可以重命名以便更清晰：
In [185]: col_mapping = {"description" : "food",
.....: "group" : "fgroup"}
In [186]: info = info.rename(columns=col_mapping, copy=False)
In [187]: info.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 6636 entries, 0 to 6635
Data columns (total 4 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 food 6636 non-null object
1 fgroup 6636 non-null object
2 id 6636 non-null int64
3 manufacturer 5195 non-null object
dtypes: int64(1), object(3)
memory usage: 207.5+ KB
In [188]: col_mapping = {"description" : "nutrient",
.....: "group" : "nutgroup"}
In [189]: nutrients = nutrients.rename(columns=col_mapping, copy=False)
In [190]: nutrients
Out[190]:
value units nutrient nutgroup id
0 25.180 g Protein Composition 1008
1 29.200 g Total lipid (fat) Composition 1008
2 3.060 g Carbohydrate, by difference Composition 1008
3 3.280 g Ash Other 1008
4 376.000 kcal Energy Energy 1008
... ... ... ... ... ...
389350 0.000 mcg Vitamin B-12, added Vitamins 43546
389351 0.000 mg Cholesterol Other 43546
389352 0.072 g Fatty acids, total saturated Other 43546
389353 0.028 g Fatty acids, total monounsaturated Other 43546
389354 0.041 g Fatty acids, total polyunsaturated Other 43546
[375176 rows x 5 columns]
完成所有这些后，我们准备将info与nutrients合并：
In [191]: ndata = pd.merge(nutrients, info, on="id")
In [192]: ndata.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 375176 entries, 0 to 375175
Data columns (total 8 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 value 375176 non-null float64
1 units 375176 non-null object
2 nutrient 375176 non-null object
3 nutgroup 375176 non-null object
4 id 375176 non-null int64
5 food 375176 non-null object
6 fgroup 375176 non-null object
7 manufacturer 293054 non-null object
dtypes: float64(1), int64(1), object(6)
memory usage: 22.9+ MB
In [193]: ndata.iloc[30000]
Out[193]:
value 0.04
units g
nutrient Glycine
nutgroup Amino Acids
id 6158
food Soup, tomato bisque, canned, condensed
fgroup Soups, Sauces, and Gravies
manufacturer
Name: 30000, dtype: object
现在我们可以制作按食物组和营养类型中位数值的图表（请参见各食物组的锌中位数值）：
In [195]: result = ndata.groupby(["nutrient", "fgroup"])["value"].quantile(0.5)
In [196]: result["Zinc, Zn"].sort_values().plot(kind="barh")
图 13.11：各食物组的锌中位数值
使用idxmax或argmax Series 方法，您可以找到每种营养素中最密集的食物。在 Jupyter 单元格中运行以下内容：
by_nutrient = ndata.groupby(["nutgroup", "nutrient"])
def get_maximum(x):
return x.loc[x.value.idxmax()]
max_foods = by_nutrient.apply(get_maximum)[["value", "food"]]
# make the food a little smaller
max_foods["food"] = max_foods["food"].str[:50]
生成的 DataFrame 太大，无法在书中显示；这里只有"Amino Acids"营养组：
In [198]: max_foods.loc["Amino Acids"]["food"]
Out[198]:
nutrient
Alanine Gelatins, dry powder, unsweetened
Arginine Seeds, sesame flour, low-fat
Aspartic acid Soy protein isolate
Cystine Seeds, cottonseed flour, low fat (glandless)
Glutamic acid Soy protein isolate
Glycine Gelatins, dry powder, unsweetened
Histidine Whale, beluga, meat, dried (Alaska Native)
Hydroxyproline KENTUCKY FRIED CHICKEN, Fried Chicken, ORIGINAL RE
Isoleucine Soy protein isolate, PROTEIN TECHNOLOGIES INTERNAT
Leucine Soy protein isolate, PROTEIN TECHNOLOGIES INTERNAT
Lysine Seal, bearded (Oogruk), meat, dried (Alaska Native
Methionine Fish, cod, Atlantic, dried and salted
Phenylalanine Soy protein isolate, PROTEIN TECHNOLOGIES INTERNAT
Proline Gelatins, dry powder, unsweetened
Serine Soy protein isolate, PROTEIN TECHNOLOGIES INTERNAT
Threonine Soy protein isolate, PROTEIN TECHNOLOGIES INTERNAT
Tryptophan Sea lion, Steller, meat with fat (Alaska Native)
Tyrosine Soy protein isolate, PROTEIN TECHNOLOGIES INTERNAT
Valine Soy protein isolate, PROTEIN TECHNOLOGIES INTERNAT
Name: food, dtype: object
13.5 2012 年联邦选举委员会数据库
美国联邦选举委员会（FEC）发布了有关政治竞选捐款的数据。这包括捐助者姓名、职业和雇主、地址以及捐款金额。2012 年美国总统选举的捐款数据作为一个 150 兆字节的 CSV 文件P00000001-ALL.csv可用（请参阅本书的数据存储库），可以使用pandas.read_csv加载：
In [199]: fec = pd.read_csv("datasets/fec/P00000001-ALL.csv", low_memory=False)
In [200]: fec.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 1001731 entries, 0 to 1001730
Data columns (total 16 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 cmte_id 1001731 non-null object
1 cand_id 1001731 non-null object
2 cand_nm 1001731 non-null object
3 contbr_nm 1001731 non-null object
4 contbr_city 1001712 non-null object
5 contbr_st 1001727 non-null object
6 contbr_zip 1001620 non-null object
7 contbr_employer 988002 non-null object
8 contbr_occupation 993301 non-null object
9 contb_receipt_amt 1001731 non-null float64
10 contb_receipt_dt 1001731 non-null object
11 receipt_desc 14166 non-null object
12 memo_cd 92482 non-null object
13 memo_text 97770 non-null object
14 form_tp 1001731 non-null object
15 file_num 1001731 non-null int64
dtypes: float64(1), int64(1), object(14)
memory usage: 122.3+ MB
注意：
有几个人要求我将数据集从 2012 年选举更新到 2016 年或 2020 年选举。不幸的是，联邦选举委员会提供的最新数据集变得更大更复杂，我决定在这里使用它们会分散我想要说明的分析技术。
数据框中的一个示例记录如下：
In [201]: fec.iloc[123456]
Out[201]:
cmte_id C00431445
cand_id P80003338
cand_nm Obama, Barack
contbr_nm ELLMAN, IRA
contbr_city TEMPE
contbr_st AZ
contbr_zip 852816719
contbr_employer ARIZONA STATE UNIVERSITY
contbr_occupation PROFESSOR
contb_receipt_amt 50.0
contb_receipt_dt 01-DEC-11
receipt_desc NaN
memo_cd NaN
memo_text NaN
form_tp SA17A
file_num 772372
Name: 123456, dtype: object
您可能会想到一些方法来开始切片和切块这些数据，以提取有关捐赠者和竞选捐款模式的信息统计。我将展示一些应用本书中技术的不同分析方法。
您会发现数据中没有政党隶属关系，因此添加这些信息会很有用。您可以使用unique获取所有唯一的政治候选人列表：
In [202]: unique_cands = fec["cand_nm"].unique()
In [203]: unique_cands
Out[203]:
array(['Bachmann, Michelle', 'Romney, Mitt', 'Obama, Barack',
"Roemer, Charles E. 'Buddy' III", 'Pawlenty, Timothy',
'Johnson, Gary Earl', 'Paul, Ron', 'Santorum, Rick',
'Cain, Herman', 'Gingrich, Newt', 'McCotter, Thaddeus G',
'Huntsman, Jon', 'Perry, Rick'], dtype=object)
In [204]: unique_cands[2]
Out[204]: 'Obama, Barack'
表示政党隶属关系的一种方法是使用字典：¹
parties = {"Bachmann, Michelle": "Republican",
"Cain, Herman": "Republican",
"Gingrich, Newt": "Republican",
"Huntsman, Jon": "Republican",
"Johnson, Gary Earl": "Republican",
"McCotter, Thaddeus G": "Republican",
"Obama, Barack": "Democrat",
"Paul, Ron": "Republican",
"Pawlenty, Timothy": "Republican",
"Perry, Rick": "Republican",
"Roemer, Charles E. 'Buddy' III": "Republican",
"Romney, Mitt": "Republican",
"Santorum, Rick": "Republican"}
现在，使用这个映射和 Series 对象上的map方法，您可以从候选人姓名计算一个政党数组：
In [206]: fec["cand_nm"][123456:123461]
Out[206]:
123456 Obama, Barack
123457 Obama, Barack
123458 Obama, Barack
123459 Obama, Barack
123460 Obama, Barack
Name: cand_nm, dtype: object
In [207]: fec["cand_nm"][123456:123461].map(parties)
Out[207]:
123456 Democrat
123457 Democrat
123458 Democrat
123459 Democrat
123460 Democrat
Name: cand_nm, dtype: object
# Add it as a column
In [208]: fec["party"] = fec["cand_nm"].map(parties)
In [209]: fec["party"].value_counts()
Out[209]:
party
Democrat 593746
Republican 407985
Name: count, dtype: int64
一些数据准备要点。首先，这些数据包括捐款和退款（负捐款金额）：
In [210]: (fec["contb_receipt_amt"] > 0).value_counts()
Out[210]:
contb_receipt_amt
True 991475
False 10256
Name: count, dtype: int64
为简化分析，我将限制数据集为正捐款：
In [211]: fec = fec[fec["contb_receipt_amt"] > 0]
由于巴拉克·奥巴马和米特·罗姆尼是主要的两位候选人，我还将准备一个只包含对他们竞选活动的捐款的子集：
In [212]: fec_mrbo = fec[fec["cand_nm"].isin(["Obama, Barack", "Romney, Mitt"])]
按职业和雇主的捐款统计
按职业捐款是另一个经常研究的统计数据。例如，律师倾向于向民主党捐款更多，而商业高管倾向于向共和党捐款更多。您没有理由相信我；您可以在数据中自己看到。首先，可以使用value_counts计算每个职业的总捐款数：
In [213]: fec["contbr_occupation"].value_counts()[:10]
Out[213]:
contbr_occupation
RETIRED 233990
INFORMATION REQUESTED 35107
ATTORNEY 34286
HOMEMAKER 29931
PHYSICIAN 23432
INFORMATION REQUESTED PER BEST EFFORTS 21138
ENGINEER 14334
TEACHER 13990
CONSULTANT 13273
PROFESSOR 12555
Name: count, dtype: int64
通过查看职业，您会注意到许多职业都指的是相同的基本工作类型，或者有几种相同事物的变体。以下代码片段演示了一种通过从一个职业映射到另一个职业来清理其中一些职业的技术；请注意使用dict.get的“技巧”，以允许没有映射的职业“通过”：
occ_mapping = {
"INFORMATION REQUESTED PER BEST EFFORTS" : "NOT PROVIDED",
"INFORMATION REQUESTED" : "NOT PROVIDED",
"INFORMATION REQUESTED (BEST EFFORTS)" : "NOT PROVIDED",
"C.E.O.": "CEO"
}
def get_occ(x):
# If no mapping provided, return x
return occ_mapping.get(x, x)
fec["contbr_occupation"] = fec["contbr_occupation"].map(get_occ)
我也会为雇主做同样的事情：
emp_mapping = {
"INFORMATION REQUESTED PER BEST EFFORTS" : "NOT PROVIDED",
"INFORMATION REQUESTED" : "NOT PROVIDED",
"SELF" : "SELF-EMPLOYED",
"SELF EMPLOYED" : "SELF-EMPLOYED",
}
def get_emp(x):
# If no mapping provided, return x
return emp_mapping.get(x, x)
fec["contbr_employer"] = fec["contbr_employer"].map(get_emp)
现在，您可以使用pivot_table按政党和职业对数据进行聚合，然后筛选出总捐款至少为 200 万美元的子集：
In [216]: by_occupation = fec.pivot_table("contb_receipt_amt",
.....: index="contbr_occupation",
.....: columns="party", aggfunc="sum")
In [217]: over_2mm = by_occupation[by_occupation.sum(axis="columns") > 2000000]
In [218]: over_2mm
Out[218]:
party Democrat Republican
contbr_occupation
ATTORNEY 11141982.97 7477194.43
CEO 2074974.79 4211040.52
CONSULTANT 2459912.71 2544725.45
ENGINEER 951525.55 1818373.70
EXECUTIVE 1355161.05 4138850.09
HOMEMAKER 4248875.80 13634275.78
INVESTOR 884133.00 2431768.92
LAWYER 3160478.87 391224.32
MANAGER 762883.22 1444532.37
NOT PROVIDED 4866973.96 20565473.01
OWNER 1001567.36 2408286.92
PHYSICIAN 3735124.94 3594320.24
PRESIDENT 1878509.95 4720923.76
PROFESSOR 2165071.08 296702.73
REAL ESTATE 528902.09 1625902.25
RETIRED 25305116.38 23561244.49
SELF-EMPLOYED 672393.40 1640252.54
这些数据以条形图形式更容易查看（"barh"表示水平条形图；请参见按职业和政党分组的总捐款）：
In [220]: over_2mm.plot(kind="barh")
图 13.12：按职业分组的政党总捐款
您可能对捐赠最多的职业或向奥巴马和罗姆尼捐款最多的公司感兴趣。为此，您可以按候选人姓名分组，并使用本章早期的top方法的变体：
def get_top_amounts(group, key, n=5):
totals = group.groupby(key)["contb_receipt_amt"].sum()
return totals.nlargest(n)
然后按职业和雇主进行汇总：
In [222]: grouped = fec_mrbo.groupby("cand_nm")
In [223]: grouped.apply(get_top_amounts, "contbr_occupation", n=7)
Out[223]:
cand_nm contbr_occupation
Obama, Barack RETIRED 25305116.38
ATTORNEY 11141982.97
INFORMATION REQUESTED 4866973.96
HOMEMAKER 4248875.80
PHYSICIAN 3735124.94
LAWYER 3160478.87
CONSULTANT 2459912.71
Romney, Mitt RETIRED 11508473.59
INFORMATION REQUESTED PER BEST EFFORTS 11396894.84
HOMEMAKER 8147446.22
ATTORNEY 5364718.82
PRESIDENT 2491244.89
EXECUTIVE 2300947.03
C.E.O. 1968386.11
Name: contb_receipt_amt, dtype: float64
In [224]: grouped.apply(get_top_amounts, "contbr_employer", n=10)
Out[224]:
cand_nm contbr_employer
Obama, Barack RETIRED 22694358.85
SELF-EMPLOYED 17080985.96
NOT EMPLOYED 8586308.70
INFORMATION REQUESTED 5053480.37
HOMEMAKER 2605408.54
SELF 1076531.20
SELF EMPLOYED 469290.00
STUDENT 318831.45
VOLUNTEER 257104.00
MICROSOFT 215585.36
Romney, Mitt INFORMATION REQUESTED PER BEST EFFORTS 12059527.24
RETIRED 11506225.71
HOMEMAKER 8147196.22
SELF-EMPLOYED 7409860.98
STUDENT 496490.94
CREDIT SUISSE 281150.00
MORGAN STANLEY 267266.00
GOLDMAN SACH & CO. 238250.00
BARCLAYS CAPITAL 162750.00
H.I.G. CAPITAL 139500.00
Name: contb_receipt_amt, dtype: float64
将捐款金额分桶
分析这些数据的一个有用方法是使用cut函数将捐助金额分成不同的桶：
In [225]: bins = np.array([0, 1, 10, 100, 1000, 10000,
.....: 100_000, 1_000_000, 10_000_000])
In [226]: labels = pd.cut(fec_mrbo["contb_receipt_amt"], bins)
In [227]: labels
Out[227]:
411 (10, 100]
412 (100, 1000]
413 (100, 1000]
414 (10, 100]
415 (10, 100]
...
701381 (10, 100]
701382 (100, 1000]
701383 (1, 10]
701384 (10, 100]
701385 (100, 1000]
Name: contb_receipt_amt, Length: 694282, dtype: category
Categories (8, interval[int64, right]): [(0, 1] < (1, 10] < (10, 100] < (100, 100
0] <
(1000, 10000] < (10000, 100000] < (10000
0, 1000000] <
(1000000, 10000000]]
然后，我们可以按姓名和 bin 标签对 Obama 和 Romney 的数据进行分组，以获得按捐款大小分组的直方图：
In [228]: grouped = fec_mrbo.groupby(["cand_nm", labels])
In [229]: grouped.size().unstack(level=0)
Out[229]:
cand_nm Obama, Barack Romney, Mitt
contb_receipt_amt
(0, 1] 493 77
(1, 10] 40070 3681
(10, 100] 372280 31853
(100, 1000] 153991 43357
(1000, 10000] 22284 26186
(10000, 100000] 2 1
(100000, 1000000] 3 0
(1000000, 10000000] 4 0
这些数据显示，奥巴马收到的小额捐款数量明显多于罗姆尼。您还可以对捐款金额进行求和，并在桶内进行归一化，以可视化每个候选人每个大小的总捐款的百分比（每个捐款大小收到的候选人总捐款的百分比显示了结果图）：
In [231]: bucket_sums = grouped["contb_receipt_amt"].sum().unstack(level=0)
In [232]: normed_sums = bucket_sums.div(bucket_sums.sum(axis="columns"),
.....: axis="index")
In [233]: normed_sums
Out[233]:
cand_nm Obama, Barack Romney, Mitt
contb_receipt_amt
(0, 1] 0.805182 0.194818
(1, 10] 0.918767 0.081233
(10, 100] 0.910769 0.089231
(100, 1000] 0.710176 0.289824
(1000, 10000] 0.447326 0.552674
(10000, 100000] 0.823120 0.176880
(100000, 1000000] 1.000000 0.000000
(1000000, 10000000] 1.000000 0.000000
In [234]: normed_sums[:-2].plot(kind="barh")
图 13.13：每个捐款大小收到的候选人总捐款的百分比
我排除了两个最大的桶，因为这些不是个人捐款。
这种分析可以以许多方式进行细化和改进。例如，您可以按捐赠人姓名和邮政编码对捐款进行汇总，以调整给出许多小额捐款与一笔或多笔大额捐款的捐赠者。我鼓励您自己探索数据集。
按州的捐款统计
我们可以通过候选人和州对数据进行汇总：
In [235]: grouped = fec_mrbo.groupby(["cand_nm", "contbr_st"])
In [236]: totals = grouped["contb_receipt_amt"].sum().unstack(level=0).fillna(0)
In [237]: totals = totals[totals.sum(axis="columns") > 100000]
In [238]: totals.head(10)
Out[238]:
cand_nm Obama, Barack Romney, Mitt
contbr_st
AK 281840.15 86204.24
AL 543123.48 527303.51
AR 359247.28 105556.00
AZ 1506476.98 1888436.23
CA 23824984.24 11237636.60
CO 2132429.49 1506714.12
CT 2068291.26 3499475.45
DC 4373538.80 1025137.50
DE 336669.14 82712.00
FL 7318178.58 8338458.81
如果您将每一行都除以总捐款金额，您将得到每位候选人每个州的总捐款相对百分比：
13.6 结论
在这本书第一版出版以来的 10 年里，Python 已经成为数据分析中流行和广泛使用的语言。您在这里所学习的编程技能将在未来很长一段时间内保持相关性。希望我们探讨过的编程工具和库能够为您提供帮助。
我们已经到达了这本书的结尾。我在附录中包含了一些您可能会发现有用的额外内容。
- 这做出了一个简化的假设，即 Gary Johnson 是共和党人，尽管后来成为了自由党候选人。
