Python学习笔记（9-2）：数据分析——数据统计

文章导读

 - 课程难度：★★★☆☆
 - 重要度：★★★☆☆
 - 预计学习时间：1.5小时
 - 简介：本节主要讲解了python数据分析中数据的统计及检验方法，包括：（1）基于numpy和pandas方法对数据进行基本的统计；（2）多种分布检验、参数检验、非参数检验、协方差与相关性的计算方法；（3）基于random包、numpy.random包和scipy.stats包的随机采样；（4）基于numpy.polynomial和scipy.stats的分布与回归拟合。
 - 重点涉及的函数和内容：mstats.normaltest()方法判断数据是否服从正态分布、mstats.chisquare()方法实现卡方检验、mstats.pearsonr()方法计算连续数据的pearson相关系数、mstats.spearmanr()方法计算离散数据的spearman相关系数、np.random.randint()、np.random.randn()。

一、基础统计

我们首先介绍numpy数组的统计方式，之后也对dataframe在统计方面的一些成员方法进行介绍。这里，由于numpy和pandas不直接支持众数的计算，我们会引入scipy下的stats包。

from scipy import stats
import pandas as pd
import numpy as np

需要准备的数据：

frame = pd.DataFrame({"a": [1, 2, 3, 4],"b": [4, 5, 6, 7],"c": [1, 2, 3, 4]} )
array = np.array(frame)
frame_nan = frame.replace(1, float('nan'))
array_nan = np.array(frame_nan)

x = np.array([100.45, 100.02, 100.0,  99.97, 100.05])
y = np.array([1.639, 0.578, 1.031, 1.031, 0.59]) 
x_nan = np.array([100.45, 100.02, np.nan,  99.97, 100.05])
y_nan = np.array([1.639, np.nan, 1.031, np.nan, 0.59])

xy_ary = np.array([x, y]) 
xy_ary_nan = np.array([x_nan, y_nan])

通常地，在处理多维（这里以二维为例）数组的时候，我们需要传递计算的轴向。为了直观，这里在部分情况下指定axis = 0（按列），在部分情况下指定axis = 1（按行）。

1、平均数、中位数、众数

求平均数、中位数、众数的方法及示例如下：

np.mean(array, axis = 0) 			# 平均数
np.median(array, axis = 0) 		    # 中位数
stats.mode(xy_ary, axis = 1) 		# 众数，按行统计，这里用xy_ary作为例子

输出结果如下：

array([2.5, 5.5, 2.5])
array([2.5, 5.5, 2.5])
ModeResult(mode=array([[99.97 ],
       [ 1.031]]), count=array([[1],
       [2]]))

对于众数的计算函数，其返回的结果包括两个部分，一是每组元素的值，二是各个元素的出现次数。因此，根据返回值继续访问所需的众数值，需要向下继续访问各种元素。因为两个部分是通过各自只包含一个元素的列表存储的，因此访问上述二维数组中，某一维度的众数需要按如下方式访问：（对于一个向量是2个0）

stats.mode(xy_ary, axis = 1)[0][0][0] 	# 另一个是stats.mode(xy_ary, axis = 1)[0][1][0]

运算结果如下：

99.97

对应于mean、median，我们都有对应的忽略缺失值的方法，即nanmean和nanmedian。对于众数已经是忽略了缺失值的了，方法如下：

np.nanmean(array_nan, axis = 0)
np.nanmedian(array_nan, axis = 0)
stats.mode(array_nan, axis = 0)

输出结果如下：

array([3. , 5.5, 3. ])
array([3. , 5.5, 3. ])
ModeResult(mode=array([[2., 4., 2.]]), count=array([[1, 1, 1]]))

2、最大值、最小值、极差、方差、标准差、分位数

（1）求最大值、最小值等各类有关数据分布的统计函数如下：

np.max(xy_ary_nan, axis = 0) 				# 最大值，等同于np.amax()
np.min(xy_ary_nan, axis = 0) 				# 最小值，等同于np.amin()
np.ptp(xy_ary_nan, axis = 0) 				# 极差，最大值减去最小值，名字特殊，
									        # 不写成range可能因为range用作序列函数
np.var(xy_ary_nan, axis = 0) 				# 方差
np.std(xy_ary_nan, axis = 0) 				# 标准差
np.percentile(array_nan, axis = 0, q = 90) 	# 分位数，在0-100中取，100为最大、0为最小
np.quantile(array_nan, axis = 0, q = 0.9) 	# 分位数，但是是在0-1中取，1为最大、0为最小

除了极差之外，上述函数同样有对应的忽略空值的方法：

np.nanmax(xy_ary_nan, axis = 0) 
np.nanmin(xy_ary_nan, axis = 0) 
np.nanstd(xy_ary_nan, axis = 0)
np.nanvar(xy_ary_nan, axis = 0)
np.nanpercentile(array_nan, axis = 0, q = 90) 
np.nanquantile(array_nan, axis = 0, q = 0.9)

最后一部分是一些补充：

np.average(xy_ary, axis = 0, weights=[1, 3]) 	
# 加权平均（平均是mean），要求轴的长度和权重向量长度等同。权重各项目的和不要求为1

np.sum(xy_ary, axis = 1) 					# 求和，按行统计
np.prod(xy_ary, axis = 1)					# 求乘积，按行统计
np.nansum(xy_ary_nan, axis = 0) 			# 忽略缺失值的求和，
np.nanprod(xy_ary_nan, axis = 0)			# 忽略缺失值的乘积

需要注意的是，我们这里numpy常用的统计里不支持峰度和偏度的计算，stats包支持偏度计算，但没有峰度计算，这个在pandas中可以实现。

3、基于numpy方法统计DataFrame

上述介绍的是numpy方法的统计，传入函数的参数也是numpy数组。实际上，我们也可以直接传递Series或DataFrame。当传入的参数是DataFrame时，对应mean、max、min、var、std五个函数，其结果是可以返回Series序列的。
这五个函数，对于参数axis，仍然是传递0为对列进行运算，传递1为对行进行运算，而且此时是已经忽略掉缺失值的计算结果，比较如下三种计算结果：：

np.mean(array_nan, axis = 0) 			
np.nanmean(array_nan, axis = 0) 
np.mean(frame_nan, axis = 0)

输出结果如下：

array([nan, 5.5, nan])
array([3. , 5.5, 3. ])
a    3.0
b    5.5
c    3.0
dtype: float64

相应地，dataframe的五个函数使用如下，在此为了使篇幅简洁，不展示运算结果：

np.mean(frame_nan, axis = 0) 		    # 平均值
np.max(frame_nan, axis = 0) 			# 最大值
np.min(frame_nan, axis = 0) 			# 最小值
np.var(frame_nan, axis = 0) 			# 方差
np.std(frame_nan, axis = 0) 			# 标准差

4、基于pandas方法进行统计

如果对dataframe进行处理，可首先考虑其内置方法，举一些基本的例子如下：

 frame_nan.count() 			    # 非零元素的计数
 frame_nan.sum() 				# 求和
 frame_nan.mean() 			    # 平均值
 frame_nan.median() 			# 中位数
 frame_nan.min() 				# 最小值
 frame_nan.max() 				# 最大值
 frame_nan.std() 				# 标准差
 frame_nan.var() 				# 方差
 frame_nan.skew() 				# 偏度
 frame_nan.kurt() 				# 峰度
 frame_nan.quantile() 			# 分位数，默认值为0.5，在0-1中间取
 frame_nan.cov()				# 协方差，返回的是dataframe
 frame_nan.corr() 				# 相关系数，返回dataframe

当然，如果只是想要查看dataframe的基本情况，我们就直接用.describe()方法即可：

frame_nan.describe()

二、统计检验及相关性

统计检验基于的是scipy.stats.mstats包，这个包包含了各种各样的统计检验方法。

from scipy.stats import mstats

我们这里依然准备一些数据：

frame_ori = pd.read_csv('lesson9\\data_test.csv').replace(0, 1)
frame = frame_ori.drop('TOOL', axis = 1)
array = np.array(frame)

x_ary1 = np.array([77, 80, 82, 90, 83, 73, 84, 79])
x_ary2 = np.array([89, 88, 87, 81, 90, 86, 91, 84])
x_ary3 = np.array([90, 85, 89, 93, 88, 87, 90, 86])

1、分布检验

在进行参数检验前，我们需要进行样本的分布检验，这是因为参数检验对于待检验参数的分布有所要求。例如，t检验要求总体分布是正态分布、总体标准差未知且样本容量小于30；方差分析要求样本相互独立、分别来自正态分布总体，且样本方差相同，等等。
在此，我们对一些分布检验方法进行介绍，包括偏度检验、峰度检验和正态分布检验，分别是计算偏度、峰度以及判断一组样本是否拥有近似正态分布的偏度与峰度。为此，我们构建一些数据：

x_ary_value = np.array(frame_ori['Value'])
xtmp1 = np.array([12] + [20] * 30) 		# 几乎所有元素一样的样本集合
xtmp2 = np.array([20] * 30) 				# 所有元素完全一样的样本集合

（1）计算峰度 mstats.kurtosis()

利用mstats.kurtosis()函数计算一组样本的峰度，样本的数目建议至少大于20。该函数有一个fisher参数，为True指示运用Fisher的峰度定义，为False时运用Pearson的定义，默认为True。
如果一组样本全部相同表示异常，Fisher定义下返回-3，Pearson定义下返回0。

mstats.kurtosis(x_ary_value, fisher = True) 
mstats.kurtosis(xtmp1, fisher = True) 
mstats.kurtosis(xtmp2, fisher = True) 

各自的运算结果如下：

-0.4234960480597709 
26.0333333333333 
-3.0

（2）计算峰度 mstats.kurtosistest()

我们利用它计算一组样本是否拥有近似正态分布的峰度：

mstats.kurtosistest(x_ary_value) 	# p > 0.05，认为不能拒绝原假设
								# 即可以认为样本拥有近似正态分布的峰度
mstats.kurtosistest(xtmp1) 		# p < 0.05，认为应拒绝原假设
								# 即不能认为样本拥有近似正态分布的峰度
mstats.kurtosistest(xtmp2) 		# 当所有元素全部一样时p值为1，并报告异常

（3）计算峰度 mstats.skew()

我们利用mstats.skew()函数计算一组样本的峰度，样本的数目建议至少大于20。如果一组样本全部相同表示异常，将会返回0。

mstats.skew(x_ary_value)
mstats.skew(xtmp1)
mstats.skew(xtmp2)

（4）计算偏度 mstats.skewtest()

与mstats.skew()相配套的是mstats.skewtest()，我们利用它计算一组样本是否拥有近似正态分布的偏度：

mstats.skewtest(x_ary_value) 		# p > 0.05，认为不能拒绝原假设，
								# 即可以认为样本拥有近似正态分布的偏度
mstats.skewtest(xtmp1) 			# p < 0.05，认为应拒绝原假设，
								# 即不能认为样本拥有近似正态分布的偏度
mstats.skewtest(xtmp2) 			# p > 0.05，认为不能拒绝原假设，
								# 即可以认为样本拥有近似正态分布的偏度
								# （意味着数据全部相等时能够通过偏度检验）

（5）判断数据是否符合正态分布 mstats.normaltest()

mstats.normaltest(x_ary_value) 		
# p > 0.05，认为不能拒绝原假设，即可以认为样本近似正态分布

mstats.normaltest(xtmp1) 			
# p < 0.05，认为不能拒绝原假设，即不能认为样本近似正态分布

mstats.normaltest(xtmp2) 
# 当所有元素全部一样时p值为1，并报告异常

输出结果如下：

NormaltestResult(statistic=0.2033283777680402, pvalue=0.903332849953003) 
NormaltestResult(statistic=74.2291626757502, pvalue=7.609254797941124e-17) 

NormaltestResult(statistic=masked, pvalue=1.0)
【异常】RuntimeWarning: divide by zero encountered in true_divide
  term2 = ma.power((1-2.0/A)/denom,1/3.0)

2、参数检验 - T检验、方差分析

（1）T检验

t检验本身包括单样本t检验、双样本t检验（二者统称独立样本t检验）和配对样本t检验。
单样本t检验：ttest_onesamp()函数
它传入两个参数a和popmean，分别是待执行t检验的序列和原假设的平均值。当检验结果大于0.05时，我们认为原假设通过，可以认为原序列a的平均值等于给出的popmean。
例如我们对x_ary1做两次t检验，备择假设分别为70和80：

mstats.ttest_onesamp(a = x_ary1, popmean = 70) 
# p小于0.05，应拒绝原假设，不能认为序列的平均值等于70

mstats.ttest_onesamp(a = x_ary1, popmean = 80) 	
# p大于0.05，不能拒绝原假设，可以认为序列的平均值等于80

mstats.ttest_1samp(a = x_ary1, popmean = 80) 		
# ttest_onesamp的同名方法

各自的运算结果如下：

Ttest_1sampResult(statistic=6.135506861249888, pvalue=0.0004742303023758493) 
Ttest_1sampResult(statistic=0.5577733510227171, pvalue=0.5943793566536572) 
Ttest_1sampResult(statistic=0.5577733510227171, pvalue=0.5943793566536572)

双样本t检验：ttest_ind()函数（ind即independent）
它传入两个参数a和b，分别是两个待比较的序列，这两个序列可以不等长。当检验结果大于0.05时，我们认为原假设通过，可以认为原序列a和b的平均值相同。
例如我们对x_ary1、x_ary2和x_ary3做三次t检验：

mstats.ttest_ind(a = x_ary1, b = x_ary2) 	
# p小于0.05，认为应拒绝原假设，即不能认为两个序列的平均值相等

mstats.ttest_ind(a = x_ary1, b = x_ary3) 	
# p小于0.05，认为应拒绝原假设，即不能认为两个序列的平均值相等

mstats.ttest_ind(a = x_ary2, b = x_ary3)	
# p大于0.05，认为不能拒绝原假设，即可以认为两个序列的平均值相等

各自的运算结果如下：

Ttest_indResult(statistic=-2.806243040080456, pvalue=0.014007122329600002)
Ttest_indResult(statistic=-3.7333702063075838, pvalue=0.0022253901123268985) 
Ttest_indResult(statistic=-1.0162612288412374, pvalue=0.326747038830542)

配对样本t检验：ttest_rel()函数（rel即related）
同样是传入两个参数a和b，分别是两个待比较的序列，这两个序列必须等长。当检验结果大于0.05时，我们认为原假设通过，可以认为原序列a和b对应的元素彼此相同。
例如我们对x_ary1、x_ary2和x_ary3做三次配对样本t检验：

mstats.ttest_rel(a = x_ary1, b = x_ary2) 	
# p小于0.05，认为应拒绝原假设，即不能认为两个序列的元素对应相等

mstats.ttest_rel(a = x_ary1, b = x_ary3) 	
# p小于0.05，认为应拒绝原假设，即不能认为两个序列的元素对应相等

mstats.ttest_rel(a = x_ary2, b = x_ary3) 	
# p大于0.05，认为不能拒绝原假设，即可以认为两个序列的元素对应相等

各自的运算结果如下：

Ttest_relResult(statistic=-2.517860727186759, pvalue=0.039935237881430635) 
Ttest_relResult(statistic=-5.400617248673217, pvalue=0.0010078164296789993) 
Ttest_relResult(statistic=-0.916515138991168, pvalue=0.38987959893526014)

（2）单因素方差分析

它的实现使用mstats.f_oneway()函数，其目标同样是比较各组序列（不一定只是两组）是否来自于同一个分布总体。
当检验结果大于0.05时，我们认为原假设通过，可以认为原有的2个或数个序列来源于同一个分布总体。但这里对于数个数列不属于同一个分布总体时，不会给出具体哪个数列：

mstats.f_oneway(x_ary1, x_ary2, x_ary3) 	
# p小于0.05，我们应拒绝原假设，不能认为两个序列来源于同一个分布总体

mstats.f_oneway(x_ary2, x_ary3) 			
# p大于0.05，我们不应拒绝原假设，可以认为两个序列来源于同一个分布总体

各自的运算结果如下：

F_onewayResult(statistic=8.76158940397351, pvalue=0.0017108832379214435)
F_onewayResult(statistic=1.0327868852459017, pvalue=0.326747038830542)

3、非参数检验

使用的非参数检验而不是参数检验的直接原因，是我们的样本无法满足参数检验的条件。幸而参数检验的方法都有非参数检验方法作为替代，这些方法不要求样本分布，但是相对的，一般而言相比参数检验拒绝原假设要难一些。

（1）卡方检验：mstats.chisquare()

卡方检验计算的是观测的各组的样本频数分布（而不是值分布）是否服从与已知的经验分布。
我们利用mstats.chisquare()实现卡方检验，该函数有三个重要的参数：f_obs，观测样本的频数；f_exp，样本的经验分布；自由度参数ddof：最终的自由度是k - 1 – ddof。其中，f_obs必须给出，f_exp的默认值，即默认的经验分布为各类相等，ddof默认值为0。
我们设置三个频数分布并利用卡方检验进行判别：

x_f1 = [10, 20, 30, 40, 50]
x_f2 = [8, 22, 29, 36, 52]
x_f3 = [20, 21, 49, 43, 50]

调用函数过程如下：

# 检验x_f1的观测样本分布是否服从x_f2的经验分布
# p大于0.05，表示可以认为该观测分布服从经验分布
mstats.chisquare(x_f1, x_f2) 

# 检验x_f1的观测样本分布是否服从x_f3的经验分布
# p小于0.05，表示不应认为该观测分布服从经验分布。	
mstats.chisquare(x_f1, x_f3) 

各自的运算结果如下：

Power_divergenceResult(statistic=1.2376684618063927, pvalue=0.8718586635835323) 

Power_divergenceResult(statistic=12.624268311975953, pvalue=0.01326536236177769)

（2）T检验的非参数方法：Mann-Whitney U test

在此首先介绍Mann-Whitney U test（曼-惠特尼U检验）方法，检验这两个总体的均值是否有显著的差别，可以接一个序列和一个值或是两个序列[ 这个方法要求两组样本的数目大于20，具体请参考：Mann-Whitney U test官方文档]。
其调用方法如下：

mstats.mannwhitneyu(x = x_ary1, y = x_ary2) 	
# p小于0.05，认为应拒绝原假设，即不能认为两个序列的元素对应相等

mstats.mannwhitneyu(x = x_ary1, y = x_ary3) 	
# p小于0.05，认为应拒绝原假设，即不能认为两个序列的元素对应相等

mstats.mannwhitneyu(x = x_ary2, y = x_ary2) 	
# p大于0.05，认为不能拒绝原假设，即可以认为两个序列的元素对应相等

各自的运算结果如下：

MannwhitneyuResult(statistic=10.0, pvalue=0.0237419551407684) 
MannwhitneyuResult(statistic=6.0, pvalue=0.007232452625263326) 
MannwhitneyuResult(statistic=32.0, pvalue=0.9578736202216409)

（3）单因素方差分析的几种非参数方法

① Kruskal-Wallis检验

它用于检验多组（两个或多个）样本之间的抽样结果代表的总体样本有无差别。它各组的元素的数目可以不同，但各组建议至少有5个值；同时需指出，其判别的结果并不能指示是哪个组的经验分布和其他的不同[ 可参考：Kruskal-Wallis检验的官方文档]

mstats.kruskalwallis(x_ary1, x_ary2, x_ary3) 	
# p小于0.05，我们应拒绝原假设，不能认为两序列源于同一分布总体
mstats.kruskalwallis(x_ary2, x_ary3)			
# p大于0.05，我们不应拒绝原假设，可以认为两序列源于同一分布总体
mstats.kruskal(x_ary2, x_ary3)  # kruskalwallis的同名方法

各自的运算结果如下：

KruskalResult(statistic=9.104376367614881, pvalue=0.010544106678558618)
KruskalResult(statistic=0.5468750000000059, pvalue=0.4595973861839395) 
KruskalResult(statistic=0.5468750000000059, pvalue=0.4595973861839395)

② Kolmogorov-Smirnov（柯尔莫戈洛夫-斯米尔诺夫）检验

该方法检验两组样本是否出自于同一个分布。它基于累计分布函数，调用方法和之前的类似，这个函数仅支持两个经验分布的比较；这个仅支持双侧检验[ 可参考：ks_twosamp方法的官方文档]。

stats.ks_twosamp(x_ary1, x_ary2) 		
# p小于0.05，认为应拒绝原假设，即应认为两分布不同

mstats.ks_twosamp(x_ary1, x_ary3) 		
# p小于0.05，认为应拒绝原假设，即应认为两分布不同

mstats.ks_twosamp(x_ary2, x_ary3)		
# p大于0.05，认为不应拒绝原假设，即可以认为两分布相同

mstats.ks_2samp(x_ary2, x_ary3) # ks_twosamp的同名方法

输出结果如下：

(0.625, 0.08786641394169106) 
(0.875, 0.004374982190571073) 
(0.25, 0.9639452436648751) 
(0.25, 0.9639452436648751)

③ Friedman test（弗雷德曼检测）检验

该方法也是检验两个抽样结果或是测试结果是否来源于同一个总体样本。要求至少传入三组变量，样本一一对应；网页文档中提示由于这是基于卡方分布的假设做出的方法，需要各组有超过10个值、且至少有6个组。

# p小于0.05，认为应拒绝原假设，即应认为三个分布不同来源于不同的分布总体
mstats.friedmanchisquare(x_ary1, x_ary2, x_ary3) 

运算结果如下：

FriedmanchisquareResult(statistic=9.75, pvalue=0.007635094218859962)

4、协方差与相关分析

我们常计算的相关性有三种，包括pearson、spearman、kendall-tau相关系数。相关性的原假设都是两个向量之间不存在相关关系，也即p小于0.05时应该认为有相关关系。同时需要注意的是，相关性的p值在样本数大于500时不一定可靠，但可能仍有一定解释性。
在此，重新定义之前的几个序列，为每个样本后面各追加四个值：

x_ary1_12 = np.array([77, 80, 82, 90, 83, 73, 84, 79, 82, 88, 90, 72])
x_ary2_12 = np.array([89, 88, 87, 81, 90, 86, 91, 84, 83, 90, 91, 79]) 
x_ary3_12 = np.array([90, 85, 89, 93, 88, 87, 90, 86, 81, 93, 94, 82])

（1）pearson相关系数

pearson相关系数，第一个值为相关系数，masked_array里的data为p值；与scipy.stats.pearsonr计算结果等同。要求两个变量符合正态分布，检测的是两变量之间的线性相关关系，多用于检测连续数据的相关[ 可参考：pearsonr方法的官方文档]：

mstats.pearsonr(x_ary1, x_ary2)	 	
# p > 0.05，支持原假设，认为两者之间不存在线性相关关系
mstats.pearsonr(x_ary1, x_ary3) 		
# p < 0.05，拒绝原假设，认为两者之间存在线性相关关系

各自的运算结果如下：

(-0.26504440762199816, masked_array(data=0.52582498,
             mask=False,
       fill_value=1e+20)) 
(0.6483907988391906, masked_array(data=0.08203027,
             mask=False,
       fill_value=1e+20))

（2）spearman相关系数

spearman相关系数不要求变量符合正态分布，常在不满足pearson相关性的样本分布条件时使用，也在离散（如定序）数据中使用

mstats.spearmanr(x_ary1, x_ary2) 	
# p > 0.05，支持原假设，认为两者之间不存在相关关系
mstats.spearmanr(x_ary1, x_ary3) 	
# p < 0.05，拒绝原假设，认为两者之间存在相关关系

各自的运算结果如下：

SpearmanrResult(correlation=0.11904761904761907, pvalue=masked_array(data=0.77888573,
             mask=False,
       fill_value=1e+20)) 
SpearmanrResult(correlation=0.5389318178609666, pvalue=masked_array(data=0.16811804,
             mask=False,
       fill_value=1e+20))

（3）kendalltau方法

kendalltau也是一种非参数方法，适用于两个分类变量均为有序分类的情况，用于反映分类变量的相关性[ 可参考：kendalltau方法的官方文档]

mstats.kendalltau(x_ary1, x_ary2)		 	
# p > 0.05，支持原假设，认为两者之间不存在相关关系
mstats.kendalltau(x_ary1, x_ary3) 		
# p < 0.05，拒绝原假设，认为两者之间存在相关关系

各自的运算结果如下：

KendalltauResult(correlation=0.14285714285714285, pvalue=0.6206907170753547) 
KendalltauResult(correlation=0.40006613209931935, pvalue=0.1702399546290051)

计算协方差和相关系数，就用我们的pandas和numpy的方法即可。
只是在这两个库下如果计算相关系数，无法得到p值，因此相关系数的计算还是建议采用上面的方法。在pandas.dataframe计算协方差如下，计算的协方差和相关系数是计算列与列之间的，在此请注意：

# 将此前的三列做成一个DataFrame
frame_tmp = pd.DataFrame({'x1':x_ary1_12, 'x2':x_ary2_12, 'x3':x_ary3_12})
frame_tmp.cov() # 协方差
frame_tmp.corr() # 相关系数

在numpy中计算协方差如下，可以传递一个包含多个待计算相关性的序列地列表，也可以传递两个序列，返回相关性数组。需要注意的是在二维numpy数组中，是以行为样本计算行与行之间的相关系数的：

np.cov([x_ary1_12, x_ary2_12, x_ary3_12]) # 传递一个向量列表
np.corrcoef([x_ary1_12, x_ary2_12, x_ary3_12])		

np.cov(x_ary1_12, x_ary2_12)			# 传递两个向量
np.corrcoef(x_ary1_12, x_ary2_12)

三、随机采样

1、基于random包的采样

首先我们用内置的random包，它的用法也很简明。我们先将这个包引入进来，并建立一些向量供后面随机采样时使用

import random
x_random = [0.1, 1.2, 2.3, 3.4, 4.5, 5.6, 6.7, 7.8, 8.9, 10.0]
x_randomary = np.array([0.1, 1.2, 2.3, 3.4, 4.5, 5.6, 6.7, 7.8, 8.9, 10.0])

random包自一个固定范围生成随机数的函数如下：

random.random() # 生成一个范围在[0, 1)的随机浮点数
random.uniform(10, 20) # 生成一个在[a, b)范围内的随机浮点数
random.randint(10, 20) # 生成一个在[a, b)范围内的随机整数

输出结果如下：

0.8641693146793187 
17.848882888626747 
13

random包从某个序列进行抽样的函数如下：

random.choice(x_random) 		# 从一个序列中抽取一个值
random.randrange(10, 20, 2) 	
# 从一个range序列中抽取一个值，上式等同于random.choice(range(10, 20, 2)) 

输出结果如下：

3.4 
10

接着，random打乱元素顺序的函数如下：

random.shuffle(x_random) 		
# 将原有序列打乱，没有返回值，直接在原值上修改
random.sample(x_random, 5) 		
# 将原有序列打乱后返回前n个值，不改变原值， 并不是连续采样五次

输出结果如下：

None
[7.8, 1.2, 3.4, 2.3, 8.9]

最后，我们能设定随机数的种子，当种子固定时，执行随机数生成时产出的结果就是固定的：

random.seed(100)
random.random()

输出结果如下：

0.1456692551041303

这个基本上已经能够满足我们大部分的需求，但是这个有一个典型的问题，就是它的采样全都是均匀采样。也即我们想基于某个分布进行采样是做不到的，在这种情况下我们可以运用scipy的stats包。

2、基于numpy.random包的采样

numpy.random包在简单的随机抽样下，它支持的功能与random大体是相似的。相比较，这个少一个shuffle功能，但是支持随机初始化某个形状的数组。其调用的方式大体是相同的。
随机采样函数：

np.random.rand(2, 3)				# 随机在[0, 1)范围内采样2行3列的数据
np.random.random((3, 4)) 			# 同上，随机在[0, 1)范围内采样，
								# 但这里传入的是元组
								# 另有等同的函数random.ranf和random.sample
np.random.randint(1, 10, size = (3, 4)) 	# 随机在某个区间（左闭右开）内整数采样
np.random.randn(2, 3) 				# 随机在标准正态分布内采样2行3列的数据

输出结果如下：

[[0.28948824 0.15065403 0.0682226 ]
 [0.28338021 0.77505492 0.58179074]] 

[[0.17280124 0.34809653 0.11356091 0.29881525]
 [0.51431626 0.76963018 0.88170669 0.49304087]
 [0.32429763 0.23642791 0.48491104 0.15764539]] 

[[1 3 3 5]
 [6 3 5 9]
 [9 3 9 5]] 

[[-0.06551024  3.91163134  0.25028061]
 [ 0.44382662 -0.07640477 -0.46540443]]

接下来也是在已有的一维样本，例如numpy数组中的采样：

np.random.choice(x_random, size = (3, 4)) 	# 从已有一维样本序列中进行随机采样
np.random.shuffle(x_randomary) 			# 打乱原有样本顺序，直接修改原值
np.random.permutation(x_randomary) 		# 打乱原有样本顺序，产出一个新副本
np.random.permutation(10) 				# 当传入的参数是一个数值x时，
									# 打乱的是np.arange(x)后的数组

输出结果如下：

[[ 5.6  7.8 10.   2.3]
 [ 1.2  8.9 10.   3.4]
 [ 2.3  7.8  4.5  4.5]] 

None 
[ 7.8  6.7  1.2  4.5  3.4  8.9 10.   2.3  0.1  5.6] 
[8 5 6 2 9 7 4 0 3 1]

我们也可以通过numpy设定随机采样的种子，保证每次随机采样的结果固定：

np.random.seed(100) 		# 设定随机数种子为100
np.random.rand(10)			# 在当前随机数种子下随机采样10个数

输出结果如下：

[0.54340494 0.27836939 0.42451759 0.84477613 0.00471886 0.12156912
 0.67074908 0.82585276 0.13670659 0.57509333]

最后，numpy支持几十种分布的随机采样，相对而言是很容易理解。我们以正态分布为例子，允许传入loc即平均值，scale即标准差，以及生成数组的形状元组：

np.random.normal(3, 2, (3, 2))
np.random.normal(loc = 3, scale = 2, size = (3, 2)) # 与上面的相同

输出结果如下：

[[2.62100834 3.51000289]
 [2.08394603 3.87032698]
 [1.8328099  4.63369414]] 

[[4.34544161 2.79117771]
 [1.93743925 5.05946537]
 [2.12372875 0.76336351]]

3、基于scipy.stats包的采样

随机采样在scipy的stats包下显得功能繁多，支持比numpy更多分布的采样、分布计算和拟合。所有的连续随机变量都是rv_continuous的派生类的对象，而所有的离散随机变量都是 rv_discrete的派生类的对象
对于不同的分布，大体上只有名称不同，而相应的调用方式是类似的。我们列举一些通用的比较有意义的函数名称，其对应的功能如下：
（1）rvs，pdf，cdf函数分别执行抽样、当前位置的概率分布和概率累积分布[ rvs：random variates，pdf：probability density function，cdf：cumulative distribution function]；
（2）与cdf相对应的**sf和ppf函数**[ sf：survival function，ppf：percent point function]，分别是1-cdf和cdf的反函数。同时，pdf，cdf，sf各有自己的对数形式；
（3）fit和expect函数，实现分布拟合和期望值。
导包：

from scipy import stats
import numpy as np
import pandas as pd

在正态分布下我们进行抽样使用loc和scale参数。其中，loc是我们的平均数，scale是方差，二者的默认是分别是0和1，即默认值是标准正态分布。
用rvs()函数进行抽样，这里的size是指抽样的样本个数，这里的抽样结果是基于所定义分布的抽样结果。根据这个调用方式，我们从正态分布中抽出一个样本用于后面的操作：

stats.norm.rvs(size = 10)
rv = stats.norm.rvs(loc = 3, scale = 2, size = 10) # 平均值为3，标准差为2
rv

输出结果如下：

[ 3.44479922  0.11356601  1.48729539  4.63290802  4.50088952  2.08810615  5.37924454 -0.38123365  0.2872019   0.53513097]

我们用pdf计算当前位置的概率分布，当然注意我们的loc和scale也要给出，不然它会在标准正态分布下计算概率值。同时，我们用cdf()函数计算当前位置的累计概率，用sf()函数计算当前位置的生存函数，这里的就像之前说的，cd = 1-s

stats.norm.pdf(rv, loc = 3, scale = 2)
stats.norm.cdf(rv, loc = 3, scale = 2)
stats.norm.sf(rv, loc = 3, scale = 2)

输出结果如下：

[0.19459856 0.07040272 0.14985005 0.14293229 0.15051848 0.17977854
 0.09830442 0.04777856 0.07950027 0.09333752] 

[0.58799859 0.07447968 0.224719   0.79287974 0.77350656 0.32421407
 0.88290256 0.045455   0.08748612 0.10889341]

[0.41200141 0.92552032 0.775281   0.20712026 0.22649344 0.67578593
 0.11709744 0.954545   0.91251388 0.89110659]

最后，我们用ppf将累积概率推算回原值，ppf是cdf的反运算，因此计算结果就是之前抽样的结果rv：

cd = stats.norm.cdf(rv, loc = 3, scale = 2)
pp = stats.norm.ppf(cd, loc = 3, scale = 2)
pp

输出结果如下：

[ 3.44479922  0.11356601  1.48729539  4.63290802  4.50088952  2.08810615
  5.37924454 -0.38123365  0.2872019   0.53513097]

四、分布与回归拟合

1、基于numpy.polynomial的多项式回归拟合

多项式序列是形如𝑝(𝑥) = 1 + 2𝑥 + 3𝑥^2，而Chebyshev series形如𝑝(𝑥) = 1𝑇0(𝑥) + 2𝑇1(𝑥) + 3𝑇2(𝑥)。我们这里介绍的是基于多项式序列的回归拟合工具numpy.polynomial.polynomial，对于其他序列可类比。
导包：

import numpy.polynomial.polynomial as Poly

（1）多项式对象的定义与运算

定义
可以用Polynomial函数来定义多项式对象，这个函数需要给出从0次到高次的系数序列。定义方式与输出结果如下：

Poly.Polynomial([1, 2, 3])

输出结果如下：

poly([1. 2. 3.])

另一种方式是给出p(x) = 0的根列表r，由系统基于式子(x - r[0])*(x - r[1])*...*(x - r[n-1])自行推断原多项式。调用方式与输出结果如下：

Poly.Polynomial.fromroots([0, 0, 1, 2]) 		
# 两个0代表有四个根，其中两个是0，这和写[0, 1, 2]即三个根含义是不一样的

输出结果如下：

poly([ 0.  0.  2. -3.  1.])

我们用第一种方式分别定义两个多项式，分别是1 + 2x + 3x^2和-1 + x：

p1 = Poly.Polynomial([1, 2, 3])
p2 = Poly.Polynomial([-1, 1])

对于定义的这个多项式对象，可以查看它的几个属性：

p1.coef			# 各个阶的系数
p1.domain		# 域，可以不做关注
p1.window 		# 窗口，可以不做关注

输出结果如下：

(array([ 1.,  2.,  3.]), array([-1,  1]), array([-1,  1]))

运算
这些对象之间能够进行常用的一些数学运算，包括加、减、乘、除、乘方计算的方式就和我们多项式之间的计算一样：

p1 + p2 			# (1 + 2x + 3x^2) + (-1 + x)
p1 - p2 			# (1 + 2x + 3x^2) - (-1 + x)
p1 * p2 			# (1 + 2x + 3x^2) * (-1 + x)
p1 // p2 			# (1 + 2x + 3x^2) // (-1 + x)，多项式的相除仅支持整数除法
p1 ** 2 			# (1 + 2x + 3x^2) * (1 + 2x + 3x^2)

输出结果如下：

poly([0. 3. 3.]) 
poly([2. 1. 3.]) 
poly([-1. -1. -1.  3.]) 
poly([5. 3.]) 
poly([ 1.  4. 10. 12.  9.])

多项式对象也能与数字和列表进行加减乘数运算，这种相对而言比较便利。与列表的运算等同于与基于列表建立的同类型多项式对象进行计算，数字则看做一元列表：

p1 * [-1, 1] 		# 等同于p1 * Poly.Polynomial([-1, 1])
 p1 // [-1, 1] 		# 等同于p1 // Poly.Polynomial([-1, 1])
 p1 * 2 			# 等同于p1 * [2]

输出结果如下：

poly([-1. -1. -1.  3.]) 
poly([5. 3.]) 
poly([2. 4. 6.])

（2）多项式对象的内置功能

求根.roots()
.roots()函数会返回复数的计算结果：

p1 = Poly.Polynomial([1, 2, 3])
p2 = Poly.Polynomial([-1, 1])
p1.roots()
p2.roots()

输出结果如下：

[-0.33333333-0.47140452j -0.33333333+0.47140452j]
[1.]

求积分和微分：.integ()和.deriv()

p3 = Poly.Polynomial([2, 6]) 			# (2 + 6x)
p4 = Poly.Polynomial([1, 2, 3])
p3.integ() 						# 计算积分，默认计算一次
p3.integ(2) 						# 连续计算两次积分
p4.deriv() 						# 计算微分，默认计算一次
p4.deriv(2)						# 连续计算两次微分

输出结果如下：

poly([0. 2. 3.]) 
poly([0. 0. 1. 1.]) 
poly([2. 6.]) 
poly([6.])

.polyint()和polyder()
.polyint()和polyder()，分别与.integ()和. deriv()方法相似。只是，二者不是多项式对象的成员方法，而是一个类方法，接收多项式系数列表：

Poly.polyint(p3.coef)		
Poly.polyint(p3.coef, 2)
Poly.polyder(p4.coef)
Poly.polyder(p4.coef, 2)

输出结果如下：

[0. 2. 3.] 
[0. 0. 1. 1.] 
[2. 6.] 
[6.]

（3）利用多项式对象进行拟合

拟合
我们定义两组数据，从图像上近似是一个二次曲线。

x = [37, 37.5, 38, 38.5, 39, 39.5, 40, 40.5, 41, 41.5, 42, 42.5, 43]
y = [3.4, 3, 3, 2.27, 2.1, 1.83, 1.53, 1.7, 1.8, 1.9, 2.35, 2.54, 2.9]

实际上，执行回归拟合我们用.fit()函数即可实现，拟合过程以最小二乘的形式进行拟合，拟合结果是一个多项式对象。
除了进行回归的x和y序列之外，还必须要传递一个deg参数，如果deg参数为整数a，那么指示在最高次项不大于a（包含）的范围内进行拟合。deg参数如果为一个列表，列表中的整数是我们拟合的多项式项，例如我们传递[0, 2, 3]代表用常数、x^2 和 x^3这三个项进行拟合，而对其他次项的系数固定为0。

fit1 = Poly.Polynomial.fit(x, y, deg = [1, 2, 5]) 	# 用x、x^2和x^5进行拟合
fit2 = Poly.Polynomial.fit(x, y, deg = 2) 		# 用最高次项不高于2的多项式进行拟合
fit1
fit2

输出结果如下：

poly([ 0.         -0.4329234   4.01585934  0.      0.     0.19633672]) 
poly([ 1.75132867 -0.3210989   1.49394605])

我们需要用.convert()函数将这种形式描述的多项式对象转换成原始的多项式形式。接着，可以对转换的结果利用.coef访问具体值：

fit1.convert().coef
fit2.convert().coef

输出结果如下，这就是我们所需要的可以应用在后续预测的多项式各阶系数：

[-8.20164282e+04  1.03061755e+04 -5.16654612e+02  1.29275205e+01
 -1.61594006e-01  8.07970030e-04]
 
[ 2.71623057e+02 -1.33865534e+01  1.65994006e-01]

而类似于上述函数还有一个父类的类方法，.polyfit()函数。但这个函数拟合后，只返回一个用于描述多项式的参数列表，而不返回一个多项式对象：

Poly.polyfit(x, y, deg = 2)

输出结果如下（系数分别为0次、1次、2次），和上面fit2的结果相同：

[ 2.71623057e+02 -1.33865534e+01  1.65994006e-01]

回归预测
拟合产出拟合系数之后，我们也许需要将结果用于其他值的回归预测，这时我们使用.polyval()函数，它需要接收一个带预测的自变量值（列表）以及一个系数列表：

coef_tmp = Poly.polyfit(x, y, deg = 2)	 # 获取参数列表
Poly.polyval([36, 36.5, 43.5, 44], coef_tmp)

输出结果如下：

[4.83536464 4.15937063 3.41013986 3.9791009 ]

2、基于scipy.stats的分布拟合

分布拟合用scipy.stats包，我们首先生成一个正态分布的随机样本用于拟合：

rv = stats.norm.rvs(loc = 3, scale = 2, size = 1000)

分布拟合使用fit，我们现有的随机抽样的样本结果进行正态分布拟合，返回最好的拟合参数loc和scale，输出结果如下：

stats.norm.fit(rv)

(3.0165382504405627, 2.0135950797909636)

此时拟合的结果就已经和原分布很近似了，我们进而可以拿着fit出来的参数，来预测某个新样本的概率和累积概率。
我们补充一个expect()函数，这个函数是计算对于一个操作function，计算当传入样本满足当前刻画的分布的时候，function的期望值。我们定义三个函数，紧接着用expect计算在当前样本分布下各函数的期望值：

def fun1(x):
    return x-1
def fun2(x):
    return x*2
def fun3(x):
    return x**2

当样本满足均值为3、标准差为2时，函数的期望值：

 stats.norm.expect(fun1, loc = 3, scale = 2)
 stats.norm.expect(fun2, loc = 3, scale = 2)
 stats.norm.expect(fun3, loc = 3, scale = 2)

输出结果如下：

2.0000000000000058 
6.000000000000012 
13.0

本章讲解了我们在数据分析领域常用的各种功能。下章我们会对机器学习的工具进行一个介绍，例如我们的分类、聚类、回归等，最后一章我们会对画图进行集中讲解。

---

对于缺乏Python基础的同仁，可以通过免费专栏🔥《Python学习笔记（基础）》从零开始学习Python

结语

请始终相信“Done is better than perfect” ，不要过分追求完美，即刻行动就是最好的开始， 一个模块一个模块地积累和练习，必将有所收获。
还有其他的问题欢迎在评论区留言📝！

---

[版权申明] 非商业目的注明出处可自由转载，转载请标明出处！！！
博客：butterfly_701c