3.1. Statistics in Python

作者GaëlVaroquaux

要求

要安装Python和这些依赖项,我们建议您下载Anaconda PythonEnthought Canopy,或者最好使用软件包管理器,如果您在Ubuntu或其他linux下。

也可以看看

Python中的贝叶斯统计

本章不包括贝叶斯统计工具。贝叶斯建模特别感兴趣的是PyMC,它在Python中实现了概率编程语言。

小费

为什么使用Python进行统计?

R是专用于统计的语言。Python是一种具有统计模块的通用语言。R具有比Python更多的统计分析功能,以及专门的语法。然而,当涉及到构建复杂的分析管道,混合统计与例如图像分析,文本挖掘或物理实验的控制,Python的丰富性是非常宝贵的资产。

小费

在本文档中,Python输入使用符号“>>>”表示。

免责声明:性别问题

本教程的一些示例是围绕性别问题选择的。原因是,在控制索赔真相的这些问题实际上对许多人很重要。

3.1.1. 数据表示和交互

3.1.1.1. 数据表格

我们考虑用于统计分析的设置是由一组不同的属性特征 t>描述的多个观察t3>。因此,数据可以被视为2D表或矩阵,其中列给出数据的不同属性,并且行给出观察值。例如,examples/brain_size.csv中包含的数据:

"";"Gender";"FSIQ";"VIQ";"PIQ";"Weight";"Height";"MRI_Count"
"1";"Female";133;132;124;"118";"64.5";816932
"2";"Male";140;150;124;".";"72.5";1001121
"3";"Male";139;123;150;"143";"73.3";1038437
"4";"Male";133;129;128;"172";"68.8";965353
"5";"Female";137;132;134;"147";"65.0";951545

3.1.1.2. pandas数据框架

小费

我们将在pandas.DataFrame中从pandas模块存储和操作此数据。它是电子表格表格的Python等效项。它不同于2D numpy数组,因为它具有命名列,可以按列包含不同数据类型的混合,并且具有精细的选择和关键机制。

3.1.1.2.1. Creating dataframes: reading data files or converting arrays

从CSV文件读取:使用上述CSV文件,可以观察到大脑大小,体重和智商(Willerman et al。1991),数据是数值和分类值的混合:

>>> import pandas
>>> data = pandas.read_csv('examples/brain_size.csv', sep=';', na_values=".")
>>> data
Unnamed: 0 Gender FSIQ VIQ PIQ Weight Height MRI_Count
0 1 Female 133 132 124 118 64.5 816932
1 2 Male 140 150 124 NaN 72.5 1001121
2 3 Male 139 123 150 143 73.3 1038437
3 4 Male 133 129 128 172 68.8 965353
4 5 Female 137 132 134 147 65.0 951545
...

警告

缺少值

CSV文件中缺少第二个人的权重。如果我们不指定缺失值(NA =不可用)标记,我们将无法进行统计分析。

从数组创建:A pandas.DataFrame也可以被视为1D'系列'的字典,例如数组或列表。如果我们有3个numpy数组:

>>> import numpy as np
>>> t = np.linspace(-6, 6, 20)
>>> sin_t = np.sin(t)
>>> cos_t = np.cos(t)

我们可以将它们作为pandas.DataFrame

>>> pandas.DataFrame({'t': t, 'sin': sin_t, 'cos': cos_t})  
cos sin t
0 0.960170 0.279415 -6.000000
1 0.609977 0.792419 -5.368421
2 0.024451 0.999701 -4.736842
3 -0.570509 0.821291 -4.105263
4 -0.945363 0.326021 -3.473684
5 -0.955488 -0.295030 -2.842105
6 -0.596979 -0.802257 -2.210526
7 -0.008151 -0.999967 -1.578947
8 0.583822 -0.811882 -0.947368
...

其他输入pandas可以从SQL,excel文件或其他格式输入数据。请参阅pandas文档

3.1.1.2.2. Manipulating data

数据是类似于R的数据帧的pandas.DataFrame

>>> data.shape    # 40 rows and 8 columns
(40, 8)
>>> data.columns # It has columns
Index([u'Unnamed: 0', u'Gender', u'FSIQ', u'VIQ', u'PIQ', u'Weight', u'Height', u'MRI_Count'], dtype='object')
>>> print(data['Gender']) # Columns can be addressed by name
0 Female
1 Male
2 Male
3 Male
4 Female
...
>>> # Simpler selector
>>> data[data['Gender'] == 'Female']['VIQ'].mean()
109.45

注意

要快速查看大型数据框架,请使用其describe方法:pandas.DataFrame.describe()

groupby:对分类变量的值拆分数据框:

>>> groupby_gender = data.groupby('Gender')
>>> for gender, value in groupby_gender['VIQ']:
... print((gender, value.mean()))
('Female', 109.45)
('Male', 115.25)

groupby_gender是一个功能强大的对象,在结果数据框组中暴露许多操作:

>>> groupby_gender.mean()
Unnamed: 0 FSIQ VIQ PIQ Weight Height MRI_Count
Gender
Female 19.65 111.9 109.45 110.45 137.200000 65.765000 862654.6
Male 21.35 115.0 115.25 111.60 166.444444 71.431579 954855.4

小费

使用groupby_gender上的制表符补全可查找更多。其他常见的分组函数是中值,计数(用于检查以查看不同子集中的缺失值的数量)或sum。

../../_images/sphx_glr_plot_pandas_001.png

练习

  • 整个人口的VIQ的平均值是多少?

  • 本研究中包括多少男性/女性?

    提示使用“制表符完成”来找出可以调用的方法,而不是上面示例中的“平均值”。

  • 以计数单位表示的男性和女性的MRI计数的平均值是多少?

注意

groupby_gender.boxplot用于上图(参见此示例)。

3.1.1.2.3. Plotting data

Pandas带有一些绘图工具(pandas.tools.plotting,使用场景后面的matplotlib)显示数据帧中数据的统计信息:

散射矩阵

>>> from pandas.tools import plotting
>>> plotting.scatter_matrix(data[['Weight', 'Height', 'MRI_Count']])
../../_images/sphx_glr_plot_pandas_002.png
>>> plotting.scatter_matrix(data[['PIQ', 'VIQ', 'FSIQ']])   
../../_images/sphx_glr_plot_pandas_003.png

练习

仅绘制男性的散点矩阵,仅绘制女性散点矩阵。你认为2个子群体对应于性别吗?

3.1.2. 假设检验:比较两组

对于简单的统计测试,我们将使用scipyscipy.stats子模块:

>>> from scipy import stats

也可以看看

Scipy是一个大型图书馆。有关整个库的快速摘要,请参阅scipy章节。

3.1.2.1. 学生t检验:最简单的统计检验

3.1.2.1.1. 单样本t检验:测试群体平均值的值

../../_images/two_sided.png

scipy.stats.ttest_1samp()测试数据的总体均值是否可能等于给定值(技术上,如果观察值来自给定总体均值的高斯分布)。它返回T统计量p值(参见函数的帮助):

>>> stats.ttest_1samp(data['VIQ'], 0)   
(...30.088099970..., 1.32891964...e-28)

小费

对于p值为10 ^ -28,我们可以声称IQ(VIQ测量)的总体均值不为0。

3.1.2.1.2. 2样本t检验:测试群体之间的差异

我们已经看到,在男性和女性人群中的平均VIQ是不同的。为了测试这是否显着,我们用scipy.stats.ttest_ind()进行2样本t检验:

>>> female_viq = data[data['Gender'] == 'Female']['VIQ']
>>> male_viq = data[data['Gender'] == 'Male']['VIQ']
>>> stats.ttest_ind(female_viq, male_viq)
(...-0.77261617232..., 0.4445287677858...)

3.1.2.2. 成对测试:对相同个体进行重复测量

../../_images/sphx_glr_plot_paired_boxplots_001.png

PIQ,VIQ和FSIQ给出3个IQ量度。让我们测试FISQ和PIQ是否显着不同。我们可以使用2个样本测试:

>>> stats.ttest_ind(data['FSIQ'], data['PIQ'])   
(...0.46563759638..., 0.64277250...)

这种方法的问题是它忘记了观察之间的链接:FSIQ和PIQ是在同一个人上测量的。因此,由于受试者间变异性引起的方差是混杂的,并且可以使用“配对测试”或“重复测量测试”来移除:

>>> stats.ttest_rel(data['FSIQ'], data['PIQ'])   
(...1.784201940..., 0.082172638183...)
../../_images/sphx_glr_plot_paired_boxplots_002.png

这相当于一个单样本测试的区别:

>>> stats.ttest_1samp(data['FSIQ'] - data['PIQ'], 0)   
(...1.784201940..., 0.082172638...)

T检验假设高斯误差。我们可以使用Wilcoxon符号秩检验,放宽以下假设:

>>> stats.wilcoxon(data['FSIQ'], data['PIQ'])   
(274.5, 0.106594927...)

注意

在非配对情况下的相应测试是Mann-Whitney U testscipy.stats.mannwhitneyu()

练习

  • 测试男性和女性的体重之间的差异。
  • 使用非参数统计来测试男性和女性的VIQ之间的差异。

结论:我们发现数据不支持男性和女性具有不同VIQ的假设。

3.1.3. 线性模型,多因素和方差分析

3.1.3.1. “formula”指定Python中的统计模型

3.1.3.1.1. A simple linear regression

../../_images/sphx_glr_plot_regression_001.png

给定两组观察结果,xy,我们想要测试这样的假设:yt3>。换句话说:

其中e是观测噪声。我们将使用statsmodels模块:

  1. 拟合线性模型。我们将使用最简单的策略,普通最小二乘法(OLS)。
  2. 测试coef为非零。

首先,我们根据模型生成模拟数据:

>>> import numpy as np
>>> x = np.linspace(-5, 5, 20)
>>> np.random.seed(1)
>>> # normal distributed noise
>>> y = -5 + 3*x + 4 * np.random.normal(size=x.shape)
>>> # Create a data frame containing all the relevant variables
>>> data = pandas.DataFrame({'x': x, 'y': y})

然后我们指定一个OLS模型并适合:

>>> from statsmodels.formula.api import ols
>>> model = ols("y ~ x", data).fit()

我们可以检查从拟合得出的各种统计数据:

>>> print(model.summary())  
OLS Regression Results
==========================...
Dep. Variable: y R-squared: 0.804
Model: OLS Adj. R-squared: 0.794
Method: Least Squares F-statistic: 74.03
Date: ... Prob (F-statistic): 8.56e-08
Time: ... Log-Likelihood: -57.988
No. Observations: 20 AIC: 120.0
Df Residuals: 18 BIC: 122.0
Df Model: 1
==========================...
coef std err t P>|t| [95.0% Conf. Int.]
------------------------------------------...
Intercept -5.5335 1.036 -5.342 0.000 -7.710 -3.357
x 2.9369 0.341 8.604 0.000 2.220 3.654
==========================...
Omnibus: 0.100 Durbin-Watson: 2.956
Prob(Omnibus): 0.951 Jarque-Bera (JB): 0.322
Skew: -0.058 Prob(JB): 0.851
Kurtosis: 2.390 Cond. No. 3.03
==========================...

术语:

Statsmodels使用统计术语:statsmodels中的y变量​​称为“内生”,而x变量称为外生变量。这在这里中有更详细的讨论。

为了简化,y(内生的)是你想要预测的值,而x(外生的)代表你用来做预测的特征。

练习

从上面的模型中检索估计的参数。提示:使用制表符补全找到相关属性。

3.1.3.1.2. Categorical variables: comparing groups or multiple categories

让我们回头脑大小的数据:

>>> data = pandas.read_csv('examples/brain_size.csv', sep=';', na_values=".")

我们可以使用线性模型来比较男性和女性的IQ:

>>> model = ols("VIQ ~ Gender + 1", data).fit()
>>> print(model.summary())
OLS Regression Results
==========================...
Dep. Variable: VIQ R-squared: 0.015
Model: OLS Adj. R-squared: -0.010
Method: Least Squares F-statistic: 0.5969
Date: ... Prob (F-statistic): 0.445
Time: ... Log-Likelihood: -182.42
No. Observations: 40 AIC: 368.8
Df Residuals: 38 BIC: 372.2
Df Model: 1
==========================...
coef std err t P>|t| [95.0% Conf. Int.]
-----------------------------------------------------------------------...
Intercept 109.4500 5.308 20.619 0.000 98.704 120.196
Gender[T.Male] 5.8000 7.507 0.773 0.445 -9.397 20.997
==========================...
Omnibus: 26.188 Durbin-Watson: 1.709
Prob(Omnibus): 0.000 Jarque-Bera (JB): 3.703
Skew: 0.010 Prob(JB): 0.157
Kurtosis: 1.510 Cond. No. 2.62
==========================...

指定模型的提示

强制类别:“性别”会自动检测为类别变量,因此每个不同的值都会被视为不同的实体。

可以强制使用以下方式将整数列视为分类:

>>> model = ols('VIQ ~ C(Gender)', data).fit()

拦截:我们可以使用公式中的 - 1删除截距,或者使用+ 1强制使用截距。

小费

默认情况下,statsmodel将K个可能值的类别变量视为K-1'dummy'布尔变量(最后一个级别被截取项吸收)。这几乎总是一个好的默认选择 - 但是,可以为分类变量指定不同的编码(http://statsmodels.sourceforge.net/devel/contrasts.html)。

链接到不同FSIQ和PIQ之间的t检验

为了比较不同类型的IQ,我们需要创建一个“长形”表,列出IQ,其中IQ的类型由分类变量指示:

>>> data_fisq = pandas.DataFrame({'iq': data['FSIQ'], 'type': 'fsiq'})
>>> data_piq = pandas.DataFrame({'iq': data['PIQ'], 'type': 'piq'})
>>> data_long = pandas.concat((data_fisq, data_piq))
>>> print(data_long)
iq type
0 133 fsiq
1 140 fsiq
2 139 fsiq
...
31 137 piq
32 110 piq
33 86 piq
...
>>> model = ols("iq ~ type", data_long).fit()
>>> print(model.summary())
OLS Regression Results
...
==========================...
coef std err t P>|t| [95.0% Conf. Int.]
------------------------------------------...
Intercept 113.4500 3.683 30.807 0.000 106.119 120.781
type[T.piq] -2.4250 5.208 -0.466 0.643 -12.793 7.943
...

我们可以看到,对于iq类型的效果,我们检索与之前的t检验相同的t检验值和相应的p值:

>>> stats.ttest_ind(data['FSIQ'], data['PIQ'])   
(...0.46563759638..., 0.64277250...)

3.1.3.2. 多元回归:包括多个因素

../../_images/sphx_glr_plot_regression_3d_001.png

考虑一个线性模型,用2个变量xy解释变量z

这样的模型可以在3D中看作将平面拟合到(xyz)点的云。

示例:iris数据examples/iris.csv

小费

萼片和花瓣尺寸倾向于相关:更大的花更大!但是,还有物种的系统效应吗?

../../_images/sphx_glr_plot_iris_analysis_001.png
>>> data = pandas.read_csv('examples/iris.csv')
>>> model = ols('sepal_width ~ name + petal_length', data).fit()
>>> print(model.summary())
OLS Regression Results
==========================...
Dep. Variable: sepal_width R-squared: 0.478
Model: OLS Adj. R-squared: 0.468
Method: Least Squares F-statistic: 44.63
Date: ... Prob (F-statistic): 1.58e-20
Time: ... Log-Likelihood: -38.185
No. Observations: 150 AIC: 84.37
Df Residuals: 146 BIC: 96.41
Df Model: 3
==========================...
coef std err t P>|t| [95.0% Conf. Int.]
------------------------------------------...
Intercept 2.9813 0.099 29.989 0.000 2.785 3.178
name[T.versicolor] -1.4821 0.181 -8.190 0.000 -1.840 -1.124
name[T.virginica] -1.6635 0.256 -6.502 0.000 -2.169 -1.158
petal_length 0.2983 0.061 4.920 0.000 0.178 0.418
==========================...
Omnibus: 2.868 Durbin-Watson: 1.753
Prob(Omnibus): 0.238 Jarque-Bera (JB): 2.885
Skew: -0.082 Prob(JB): 0.236
Kurtosis: 3.659 Cond. No. 54.0
==========================...

3.1.3.3. 事后假设检验:方差分析(ANOVA)

在上述虹膜实例中,我们希望测试在除去萼片宽度的影响之后花瓣长度是否与花色和维吉尼亚不同。这可以被公式化为测试上面估计的线性模型中与彩色和维吉尼亚相关的系数之间的差异(它是方差分析,ANOVA)。为此,我们在估计的参数上写一个“对比度”的向量:我们要测试“name [T.versicolor] - name [T.virginica]“F-test

>>> print(model.f_test([0, 1, -1, 0]))
<F test: F=array([[ 3.24533535]]), p=[[ 0.07369059]], df_denom=146, df_num=1>

这个差别是否重要?

练习

回到脑大小+ IQ数据,测试是否除去脑大小,身高和体重的影响后,男女的VIQ是不同的。

3.1.4. 更多可视化:seaborn用于统计探索

Seaborn将简单的统计拟合与Pandas数据框架上的绘图相结合。

Let us consider a data giving wages and many other personal information on 500 individuals (Berndt, ER. The Practice of Econometrics. 1991. NY: Addison-Wesley).

小费

工资数据的完整代码加载和绘图在对应示例中找到。

>>> print data  
EDUCATION SOUTH SEX EXPERIENCE UNION WAGE AGE RACE \
0 8 0 1 21 0 0.707570 35 2
1 9 0 1 42 0 0.694605 57 3
2 12 0 0 1 0 0.824126 19 3
3 12 0 0 4 0 0.602060 22 3
...

3.1.4.1. Pairplot:散点矩阵

我们可以很容易地使用seaborn.pairplot()来显示连续变量之间的交互,以显示散点矩阵:

>>> import seaborn
>>> seaborn.pairplot(data, vars=['WAGE', 'AGE', 'EDUCATION'],
... kind='reg')
../../_images/sphx_glr_plot_wage_data_001.png

分类变量可以绘制为色调:

>>> seaborn.pairplot(data, vars=['WAGE', 'AGE', 'EDUCATION'],
... kind='reg', hue='SEX')
../../_images/sphx_glr_plot_wage_data_002.png

Look and feel和matplotlib设置

Seaborn改变了默认的matplotlib数字,以实现更“现代”,“excel-like”的外观。它在导入时。您可以使用以下方式重置默认值:

>>> from matplotlib import pyplot as plt
>>> plt.rcdefaults()

小费

要切换回海运设置或了解海运更好的造型,请参阅海运文档的相关部分。

3.1.4.2. lmplot:绘制单变量回归

可以使用seaborn.lmplot()来绘制一个变量与其他变量之间的关系(例如工资和教育)的回归:

>>> seaborn.lmplot(y='WAGE', x='EDUCATION', data=data)  
../../_images/sphx_glr_plot_wage_data_005.png

鲁棒回归

小费

考虑到,在上面的图中,似乎有一些数据点在主云之外,它们可能是离群值,不代表人口,但驱动回归。

要计算对异常值较不敏感的回归,必须使用稳健模型这在海图中使用绘图函数中的robust=True来完成,或者在statsmodel中,通过“Robust Linear Model”替换OLS的使用,statsmodels.formula.api.rlm()

3.1.5. 测试互动

../../_images/sphx_glr_plot_wage_education_gender_001.png

男性的教育比女性的教育增加更多吗?

小费

上面的图是由两种不同的拟合。我们需要制定一个单一模型来测试整个人口的斜率差异。这是通过“交互”完成的。

>>> result = sm.ols(formula='wage ~ education + gender + education * gender',
... data=data).fit()
>>> print(result.summary())
...
coef std err t P>|t| [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept 0.2998 0.072 4.173 0.000 0.159 0.441
gender[T.male] 0.2750 0.093 2.972 0.003 0.093 0.457
education 0.0415 0.005 7.647 0.000 0.031 0.052
education:gender[T.male] -0.0134 0.007 -1.919 0.056 -0.027 0.000
==========================...
...

我们可以得出结论,教育福利男性比女性多?

接听归属讯息

  • 假设检验和p值给出效应/差异的显着性
  • 公式(带有分类变量)使您能够在数据中表达丰富的链接
  • 可视化您的数据和简单模型适合问题!
  • 调节(添加可以解释所有或部分变化的因素)是改变解释的重要建模方面。