5.1 基本统计概念

5.1.1 随机数和分布

1.rand和random_sample

rand和random_sample都是均匀分布(uniform distribution)的随机数生成函数。这两个函数除了参数略微不同之外,功能都是一样的。rand函数传入多个整数参数作为生成数组的维度。rand_sample传入一个n维的元组元素作为参数。一般来讲,推荐使用random_sample。rand主要是为了照顾习惯于使用MATLAB的用户。示例代码如下:


import numpy as np
import matplotlib.pyplot as plt

data_uniform = np.random.random_sample((10000)) 

plt.hist(data_uniform,bins=30) 

plt.show()

运行结果如图5-1所示。

图 5-1

2.randn和standard_normal

randn和standard_normal是正态分布的随机数生成函数。这两个函数的参数与rand、random_sample类似。同样地,推荐使用standard_normal。示例代码如下:


import numpy as np
import matplotlib.pyplot as plt

data_normal=np.random.standard_normal((10000))

plt.hist(data_normal,bins=30)

plt.show()

运行结果如图5-2所示。

这里生成了一个大小为10的一维随机ndarray数组。

如果想从正态分布中取样,可以使用如下命令:


data= np.random.randn(10)

其中,randn的意思是从正态分布(Normal Distribution)中取样。

3.randint和random_integers

randint和random_integers是均匀分布的整数生成函数。函数需要传入三个参数,即low、high和size。这三个参数分别代表区间的最小值和最大值,以及需要生成的数组大小。其中,size是一个n维元组数据。randint和random_integers的区别在于,randint的范围不包括最大值,而random_integers包括最大值。示例代码如下:


In[22]: x1=np.random.randint(1,10,(100))
In[23]: x1.max()
Out[23]: 9

In[24]: x2=np.random.random_integers(1,10,(100))
In[25]: x2.max()
Out[25]: 10

图 5-2

4.shuffle

shuffle可以随机打乱一个数组,并且改变此数组本身的排列。示例代码如下:


In[34]: x=np.arange(10)
In[35]: x
Out[35]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In[36]: np.random.shuffle(x)
In[37]: x
Out[37]: array([5, 9, 3, 4, 6, 1, 2, 8, 0, 7])

5.Permutation

Permutation用于返回一个打乱后的数组值,但是并不会改变传入的参数数组本身。示例代码如下:


In[40]: x=np.arange(10)
In[41]: y=np.random.permutation(x)
In[42]: x
Out[42]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In[43]: y
Out[43]: array([9, 6, 4, 8, 2, 1, 0, 7, 3, 5])

除了以上所讲的函数之外,NumPy还提供了其他大量的随机数生成函数。所有的随机数生成函数都包含在numpy.random模块中。

6.二项式分布

N次伯努利试验的结果分布即为二项分布。使用binomial(1,p)即为一次二项分布。使用binomial(1,p,n)即表示生成n维的二项分布数组,也就是伯努利分布。生成分布函数图像的代码如下:


import numpy as np
binomial = np.random.binomial(100,0.5,10000)
import seaborn as sns
sns.distplot(binomial)

运行结果如图5-3所示。

图 5-3

7.贝塔分布

考察二项分布的分布函数:

如果我们要估计参数μ,可以获得的一个最大似然估计量为:

当数据量比较小的时候,有可能会出现过拟合的问题。贝叶斯统计则可以从另外一个侧面避免这类过拟合问题。

如果读者对贝叶斯统计有所了解,就会知道贝叶斯统计实质上是在解决如下问题:

其中,Pr(p|Data)称为后验概率,Pr(p)称为先验概率,Pr(Data)是数据的边缘分布函数,Pr(Data|p)称为似然函数。可以明确的一点是,后验概率与似然函数和先验概率的乘积成正比。而边缘分布只是一个正则项,其存在的意义是让后验分布成为一个满足概率定义的数值,即其值域属于[0,1],然而一般情况下,这个边缘分布函数是很难计算的。不过,先验概率的存在,使得整个问题的解决方式都发生了变化,比如,我们用经验判断某个参数符合一个与似然函数“共轭”的先验分布函数之后,如果有了新数据,那么只需要在之前的经验判断的基础上进行一个较小的更新即可。这样就可以避免频率主义学派在最大化似然函数的流程当中仅仅通过较小的数据量得出过度拟合数据的问题。

以二项分布为例,我们可以发现,需要估计的参数与其似然函数(二项分布的)中的组合运算项无关,而仅与μm(1-μ)N-m有关。可以猜想的是,假如先验概率也是一个与μa(1-μ)b有关的分布,那么后验分布则可以很方便地计算出来了,因为只需要在μm+a(1-μ)N-m+b前面乘一个与待估计参数无关的正则项,即可获得后验分布。与二项分布有关的共轭分布即Beta分布。

beta(a,b)从Beta(a,b)分布中生成一个随机数。Beta(a,b)、rsv(n)生成一个n维的数组,每个数据都是Beat(a,b)分布中的随机数。示例代码如下:


from scipy import stats
beta = stats.beta(3,4).rvs(100000)
import seaborn as sns
sns.distplot(beta)
Out[110]: <matplotlib.axes._subplots.AxesSubplot at 0x277ce292f60>

运行结果如图5-4所示。

图 5-4

5.1.2 随机数种子

计算机生成的随机数,其实并不是真正的随机数,而是使用特定算法生成的“伪随机数”。这样,我们会面对如下两个问题。

一是,如果初始状态值一样,那么按照同样的算法得到的“随机数”结果应该是一样的。这样就不能表现出“随机”的效果。

二是,很多随机实验,有时候可能会需要再现之前的实验结果,我们正好又希望,每次实验生成的随机数都是一样的。

也正是为了解决这个问题,引入了“种子”的概念。在生成随机数之前,先定义一个“种子”,如果这个种子是一样的,那么每次随机数的生成结果也都是一样的,这样我们便可以重复再现同一个随机实验。如果“种子”不一样,那么每次随机数生成结果都不一样,这样我们就可以大量使用不同的随机数进行实验。在NumPy里面,可以使用seed函数来定义种子。调用seed(),会根据系统提供的数据进行随机初始化,也就是每次得到的随机数都会不一样。调用seed(x)的时候,计算机会根据x值进行初始化,如果x值是同一个x值,那么随机得到的结果也是一样的。示例代码如下:


In[57]: np.random.seed()
In[58]: np.random.randn()
Out[58]: 0.30163810811524266
In[59]: np.random.seed()
In[60]: np.random.randn()
Out[60]: -0.6254191344102689

In[61]: np.random.seed(1)
In[62]: np.random.randn()
Out[62]: 1.6243453636632417
In[63]: np.random.seed(1)
In[64]: np.random.randn()
Out[64]: 1.6243453636632417

5.1.3 相关系数

两个随机变量之间的相关性与独立性是统计学研究的一个非常重要的话题。先来看一个例子。

假设有两个随机变量x和y,满足如下条件:

x=0.5y+ε,其中ε~N(0,0.5)

如果从“上帝”的角度来看这两个随机变量,可以发现,x与y有一定的“线性”关联性,但它们各自却有一定的随机性。统计上用协方差来描述这种“线性”关联性,公式如下:

我们先按照上述公式生成10 000个随机数,示例代码如下:


import numpy as np
y = np.random.normal(size=10000)
x = 0.5*y + 0.5*np.random.normal(size=10000)

再计算它们之间的协方差,示例代码如下:


np.cov(x,y)
Out[80]: 
array([[0.49287893, 0.49353615],
       [0.49353615, 0.98433829]])

可以看到矩阵的对角线元素代表的是两个随机变量之间的协方差。从公式也可看出,这个数值描述了两个变量之间的线性相关性。

用协方差描述相关性存在一个问题,即会受到随机变量的量纲影响。比如,如果我们想要衡量教育年限和薪资的水平,将教育年限的单位换为天,同时将薪资水平的单位从人民币元换成角,那么协方差将会完全不同。为了避免这种量纲上的区别造成的差异,统计学中引入了皮尔逊相关系数(Pearson Correlation Coefficient)。计算公式如下:

还是用上述生成的随机数来进行计算:


np.corrcoef(x,y)
Out[6]: 
array([[1.       , 0.7058169],
       [0.7058169, 1.       ]])

可以看到,斜对角线上的数值代表了两个变量的相关系数值。这个数值区间为[-1,1],且完全不受量纲影响。

理论上来讲,只要两个随机变量存在这种线性相关性,相关系数就不会为0。而且,如果两个变量之间存在非线性相关(如quardratic、exponential、logarithmic等关系),则相关系数将无法捕捉到这种线性相关。然而可以从公式中推测的是,只要两个随机变量存在的函数关系能够被线性逼近,则两个随机变量的相关系数就不会为0。

5.1.4 基本统计量

下面就来介绍一些常用的统计量,所谓统计量,就是利用数据的函数变化,从某种维度来反映全体数据集的特征的一种函数。

1.mean

mean函数可用于计算数组的平均值。函数的第二个参数可以用来选择不同的维度进行计算。mean既可以作为一个函数,也可以作为数组的方法来调用。示例代码如下:


In[80]: x=arange(20)
In[81]: x.mean()
Out[81]: 9.5
In[82]: mean(x)
Out[82]: 9.5
In[83]: x=reshape(arange(20),(5,4))
In[84]: mean(x,0)
Out[84]: array([  8.,   9.,  10.,  11.])
In[85]: x.mean(1)
Out[85]: array([  1.5,   5.5,   9.5,  13.5,  17.5])

2.median

median函数可用于计算数组的中位数。函数的第二个参数可以用来选择不同的维度进行计算。median既可以作为一个函数,也可以作为数组的方法来调用。示例代码如下:


In[93]: x=random.randn(4,5)
In[94]: x
Out[94]: 
array([[-0.87785842,  0.04221375,  0.58281521, -1.10061918,  1.14472371],
       [ 0.90159072,  0.50249434,  0.90085595, -0.68372786, -0.12289023],
       [-0.93576943, -0.26788808,  0.53035547, -0.69166075, -0.39675353],
       [-0.6871727 , -0.84520564, -0.67124613, -0.0126646 , -1.11731035]])
In[95]: median(x)
Out[95]: -0.33232080324099667
In[96]: median(x,0)
Out[96]: array([-0.78251556, -0.11283717,  0.55658534, -0.68769431, -0.25982188])
In[97]: median(x,1)
Out[97]: array([ 0.04221375,  0.50249434, -0.39675353, -0.6871727 ])

当计算的数组大小是偶数时,median函数将返回中间两数的均值。

3.std

std函数用于计算数组的标准差。函数的第二个参数可以用来选择不同的维度进行计算。std既可以作为一个函数,也可以作为数组的方法来调用。示例代码如下:


x = np.arange(1,100,1)
x.std()
Out[12]: 28.577380332470412

4.var

var函数用于计算数组的方差。函数的第二个参数可以用来选择不同的维度进行计算。var既可以作为一个函数,也可以作为数组的方法来调用。示例代码如下:


x = np.arange(1,100,1)
x.var()
Out[13]: 816.6666666666666

5.1.5 频率分布直方图

频率分布直方图其实是对一个变量的分布密度(分布)函数进行近似估计的一个手段。考察概率密度函数与直方图的x、y轴可以得知,密度(分布)函数图像的x轴代表的是随机变量的取值,而对于离散随机变量而言,y轴代表的是对应x取值出现的概率;对于连续随机变量而言,y轴代表的是对应x取值针对其他取值的相对可能性。直方图x轴代表的是样本的取值,而y轴则代表了x在区间取值的频率,两者之间的关系就可以体现出频率分布的状况。

1.histogram

np.histogram可以用来计算一位数组的直方图数据。可以用可选参数k来定义直方图的箱体数。如果k省略不写,则默认k=10。histogram返回两个值,第一个值是k维的向量,包含了每个箱体中的样本数量;第二个值是k+1个标识箱体的端点值。示例代码如下:


import numpy as np
x = np.random.normal(size = 100)
his = np.histogram(x,bins = 10)

Out[22]: 
(array([ 4,  6, 10, 11, 14, 16, 15, 13,  7,  4], dtype=int64),
 array([-2.11227063, -1.72285629, -1.33344194, -0.9440276 , -0.55461326,
        -0.16519892,  0.22421543,  0.61362977,  1.00304411,  1.39245846,
         1.7818728 ]))

如果我们直接调用matplotlib来画图,则可以得出如图5-5所示的可视化图像,示例代码如下:


import matplotlib.pyplot as plt
plt.hist(x,bins = 10)

运行结果如图5-5所示。

图 5-5

可以看出,上述二维数组当中,第一维数组是在描述histogram的纵轴,即每个bin对应的数据个数,而第二维数组是在描述横轴,代表随机变量的取值。

2.histogram2d

histogram2d(x,y)可用于计算2维的直方图数据。可以用可选参数bins来定义直方图的箱体数。bins既可以是一个整数,也可以是一个包含两个元素的列表,分别表示各维度的箱体数。

下面来看个相关的示例,同样,首先生成两个随机数组,代码如下:


import numpy as np
x = np.random.normal(size = 100)
y = np.random.normal(size = 100)

同时调用相应函数算出histogram的数组:


(array([[0., 0., 0., 1., 0., 0., 0., 1., 0., 0.],
        [0., 0., 0., 0., 0., 1., 0., 0., 1., 0.],
        [0., 0., 2., 0., 1., 1., 2., 0., 1., 0.],
        [0., 0., 0., 1., 1., 3., 2., 1., 0., 1.],
        [1., 0., 0., 3., 4., 5., 1., 4., 0., 2.],
        [0., 1., 3., 2., 2., 3., 4., 2., 1., 0.],
        [0., 0., 1., 0., 2., 2., 7., 1., 1., 0.],
        [1., 0., 0., 6., 4., 6., 1., 0., 0., 0.],
        [0., 0., 1., 0., 2., 0., 1., 1., 0., 0.],
        [0., 0., 0., 2., 2., 0., 0., 1., 0., 0.]]),
 array([-2.84879607, -2.34851871, -1.84824134, -1.34796398, -0.84768661,
        -0.34740925,  0.15286812,  0.65314549,  1.15342285,  1.65370022,
         2.15397758]),
 array([-2.84250735, -2.29188082, -1.74125429, -1.19062776, -0.64000123,
        -0.0893747 ,  0.46125183,  1.01187836,  1.56250489,  2.11313142,
         2.66375795]))

画出相应的图,如图5-6所示。

图 5-6

如何解读画出来的图片(图5-6)?首先我们看一下生成的histogram三维数组,数组的第一维其实是一个矩阵,代表了立体图对应的x、y坐标对应的第三维坐标z值的大小,对应到图像上面便可以想象了:颜色的深浅代表了该处z值的大小,也就代表了x、y在该点附近的分布密集程度。

我们可以将数据量扩大来看一下画出来的图像可能是什么样的。示例代码如下:


import numpy as np
x = np.random.normal(size = 100000)
y = np.random.normal(size = 100000)

plt.hist2d(x,y,bins = [100,100])

画出相应的图,如图5-7所示。

图 5-7