统计-stats

时间 2019-11-11

标签统计 stats 繁體版

原文原文链接

统计-stats

SciPy的stats模块包含了多种几率分布的随机变量[1]，随机变量分为连续和离散两种。全部的连续随机变量都是rv_continuous的派生类的对象，而全部的离散随机变量都是rv_discrete的派生类的对象。html

Footnotes数组

[1]	本节中的随机变量是指几率论中的概念，不是Python中的变量

连续和离散几率分布

可使用下面的语句得到stats模块中全部的连续随机变量：dom

>>> from scipy import stats
>>> [k for k,v in stats.__dict__.items() if isinstance(v, stats.rv_continuous)]
['genhalflogistic','triang','rayleigh','betaprime', ...]

连续随机变量对象都有以下方法：函数

rvs：对随机变量进行随机取值，能够经过size参数指定输出的数组的大小。
pdf：随机变量的几率密度函数。
cdf：随机变量的累积分布函数，它是几率密度函数的积分。
sf：随机变量的生存函数，它的值是1-cdf(t)。
ppf：累积分布函数的反函数。
stat：计算随机变量的指望值和方差。
fit：对一组随机取样进行拟合，找出最适合取样数据的几率密度函数的系数。

03-scipy/scipy_stats.pyspa

几率密度函数、直方图统计和累积分布函数code

下面以正规分布为例，简单地介绍随机变量的用法。下面的语句得到缺省正规分布的随机变量的指望值和方差，咱们看到缺省状况下它是一个均值为0，方差为1的随机变量：orm

>>> stats.norm.stats()
(array(0.0), array(1.0))

norm能够像函数同样调用，经过loc和scale参数能够指定随机变量的偏移和缩放参数。对于正态分布的随机变量来讲，这两个参数至关于指定其指望值和标准差[2]：htm

>>> X = stats.norm(loc=1.0, scale=2.0)
>>> x.stats()
(array(1.0), array(4.0))

下面调用随机变量X的rvs()方法，获得包含一万次随机取样值的数组x，而后调用NumPy的mean()和var()计算此数组的均值和方差，其结果符合随机变量X的特性：对象

>>> x = X.rvs(size=10000) # 对随机变量取10000个值
>>> np.mean(x) # 指望值
1.0181259658732724
>>> np.var(x)  # 方差
4.00188661640059

咱们也可使用fit()方法对随机取样序列x进行拟合，它返回的是与随机取样值最吻合的随机变量的参数：排序

>>> stats.norm.fit(x) # 获得随机序列指望值和标准差
array([ 1.01810091,  2.00046946])

接下来比较随机变量X的几率密度函数和对数组x进行直方图统计的结果：

>>> t = np.arange(-10, 10, 0.01)
>>> pl.plot(t, X.pdf(t)) # 绘制几率密度函数的理论值
>>> p, t2 = np.histogram(x, bins=100, normed=True)
>>> t2 = (t2[:-1] + t2[1:])/2
>>> pl.plot(t2, p) # 绘制统计所获得的几率密度

其中histogram()对数组x进行直方图统计，它将数组x的取值范围分为100个区间，并统计x中的每一个值落入各个区间中的次数。histogram()返回两个数组p和t2，其中p表示各个区间的取样值出现的频数，因为normed参数为True，所以p的值是正规化以后的结果。t2表示区间，因为其中包括区间起点和终点，所以t2的长度为101。【图:正规分布的几率密度函数(左)和累积分布函数(右)】(左)显示了几率密度函数和直方图统计的结果，能够看出两者是一致的。

下面的程序绘制随机变量X的累积分布函数和数组p的累加结果，其结果如【图:正规分布的几率密度函数(左)和累积分布函数(右)】(右)所示。

>>> pl.plot(t, X.cdf(t))
>>> pl.plot(t2, np.add.accumulate(p)*(t2[1]-t2[0]))

正规分布的几率密度函数(左)和累积分布函数(右)

有些随机分布除了loc和scale参数以外，还须要额外的形状参数。例如伽玛分布可用于描述等待个独立的随机事件发生所需的时间，就是伽玛分布的形状参数。下面计算形状参数为1和2时的伽玛分布的指望值和方差：

>>> stats.gamma.stats(1.0)
(array(1.0), array(1.0))
>>> stats.gamma.stats(2.0)
(array(2.0), array(2.0))

伽玛分布的尺度参数 $\theta$ 和随机事件发生的频率相关，由scale参数指定：

>>> stats.gamma.stats(2.0, scale=2)
(array(4.0), array(8.0))

根据伽玛分布的数学定义可知其指望值为 $k \theta$ ，方差为 $k \theta^2$ 。上面的程序验证了这两个公式。

当随机分布有额外的形状参数时，它所对应的rvs()、pdf()等方法都会增长额外的参数接收形状参数。例以下面的程序调用rvs()对 $k=2,\theta=2$ 的伽玛分布取4个随机值：

>>> x = stats.gamma.rvs(2, scale=2, size=4)
>>> x
array([ 2.20814048,  3.56652153,  4.30088176,  0.68262888])

接下来调用pdf()查看上面4个随机值所对应的几率密度：

>>> stats.gamma.pdf(x, 2, scale=2)
array([ 0.18301012,  0.1498734 ,  0.12519094,  0.12130919])

也能够先建立将形状参数和尺度参数固定的随机变量，而后再调用其pdf()计算几率密度：

>>> X = stats.gamma(2, scale=2)
>>> X.pdf(x)
array([ 0.18301012,  0.1498734 ,  0.12519094,  0.12130919])

当分布函数的值域为离散时咱们称之为离散几率分布。例如投掷有六个面的骰子时，只能得到1到6的整数，所以所获得的几率分布为离散的。对于离散随机分布，一般使用几率质量函数(PMF)描述其分布状况。

在stats库中全部描述离散分布的随机变量都从rv_discrete类继承。也能够直接用rv_discrete类自定义离散几率分布。例如假设有一个不均匀的骰子，它的各点出现的几率不相等。咱们能够用下面的数组x保存骰子的全部可能值，数组p保存每一个值出现的几率：

>>> x = range(1,7)
>>> p = (0.4, 0.2, 0.1, 0.1, 0.1, 0.1)

因而能够用下面的语句定义表示这个特殊骰子的随机变量，并调用其rvs()方法投掷此骰子20次，得到符合几率p的随机数：

>>> dice = stats.rv_discrete(values=(x,p))
>>> dice.rvs(size=20)
array([2, 5, 1, 2, 1, 1, 2, 4, 1, 3, 1, 1, 4, 3, 1, 1, 1, 2, 6, 4])

Footnotes

[2]	标准差是方差的算术平方根，所以标准差为2.0时，方差为4.0

二项、泊松、伽玛分布

本节用几个实例程序对几率论中的二项分布、泊松分布以及伽玛分布进行一些实验和讨论。

二项分布是最重要的离散几率分布之一。假设有一种只有两个结果的试验，其成功几率为，那么二项分布描述了进行次这样的独立试验，成功次的几率。二项分布的几率质量函数公式以下：

$f(k;n,p) = \frac{n!}{k!(n-k)!} p^k(1-p)^{n-k}$

例如能够经过二项分布的几率质量公式计算投掷5次骰子出现3次6点的几率。投掷一次骰子，点数为6的几率(即试验成功的几率)为，试验次数为。使用二项分布的几率质量函数pmf()能够很容易计算出现次6点的几率。和几率密度函数pdf()相似，pmf()的第一个参数为随机变量的取值，后面的参数为描述随机分布所需的参数。对于二项分布来讲，参数分别为和，而取值范围则为0到之间的整数。下面的程序计算为0到6所对应的几率：

>>> stats.binom.pmf(range(6), 5, 1/6.0)
array([0.401878, 0.401878, 0.160751, 0.032150, 0.003215, 0.000129])

由结果可知：出现0或1次6点的几率为40.2%，而出现3次6点的几率为3.215%。

在二项分布中，若是试验次数很大，而每次试验成功的几率很小，其乘积比较适中时，那么试验成功的次数的几率能够用泊松分布近似描述。

在泊松分布中使用 $\lambda$ 描述单位时间(或单位面积)中随机事件的平均发生率。若是咱们将二项分布中的试验次数看做单位时间中所作的试验次数，那么它和事件出现的几率的乘积就是事件的平均发生率，即 $\lambda=n p$ 。泊松分布的几率质量函数公式以下：

$f(k;\lambda)=\frac{e^{-\lambda}\lambda^k}{k!}$

下面的程序分别计算二项分布和泊松分布的几率质量函数，其结果如【图:当n足够大时二项分布和泊松分布近似相等】所示，能够看出当足够大时，两者是十分接近的。

03-scipy/scipy_binom_poisson.py

比较二项分布和泊松分布的几率质量函数

程序中事件平均发生率 $\lambda$ 恒等于10。根据二项分布的试验次数，计算每次事件出现的几率 $p=\lambda/n$ 。程序中运算部分大体以下：

>>> _lambda = 10.0
>>> k = np.arange(20)
>>> possion = stats.poisson.pmf(k, _lambda)  # 泊松分布
>>> binom100 = stats.binom.pmf(k, 100, _lambda/100)    # 二项分布 n=100
>>> binom1000 = stats.binom.pmf(k, 1000, _lambda/1000) # 二项分布 n=1000
>>> np.max(np.abs(binom100-possion)) # 计算最大偏差
0.006755311103353312
>>> np.max(np.abs(binom1000-possion)) # n为1000时，偏差较小
0.00063017540509099912

当n足够大时二项分布和泊松分布近似相等

泊松分布适合描述单位时间内随机事件发生的次数的分布状况。例如某一设施在必定时间内的使用次数，机器出现故障的次数，天然灾害发生的次数等等。

为了加深读者对泊松分布概念的理解，下面咱们使用随机数模拟泊松分布，并与其几率质量函数进行比较，其结果如【图:模拟泊松分布】所示。图中，每秒内事件的平均发生次数为10，即 $\lambda=10$ 。其中左图的观察时间为1000秒，而右图的观察时间为50000秒。能够看出观察时间越长，每秒内事件发生的次数越符合泊松分布。

03-scipy/scipy_poisson_sim.py

模拟泊松分布

因为上面的程序中包含了许多绘图方面的代码，下面咱们直接在IPython中介绍泊松分布的模拟过程。首先定义事件发生率 $\lambda$ 和观察时间:

>>> _lambda = 10
>>> time = 10000

能够用NumPy的随机数生成函数rand()，产平生均分布于0到time之间的_lambda*time个事件所发生的时刻。因为rand()产生的是0到1之间的平均分布的随机数，所以须要对其结果扩大time倍：

>>> t = np.random.rand(_lambda*time)*time

用histogram()能够统计数组t中每秒以内的事件发生的次数count：

>>> count, time_edges = np.histogram(t, bins=time, range=(0,time))
>>> count
array([10,  9,  8, ..., 11, 10, 18])

根据泊松分布的定义，count数组中的数值的分布状况应该符合泊松分布。下面统计事件次数在0到20区间内的几率分布。当histogram()的normed参数为True而且每一个统计区间的长度为1时，其结果和几率质量函数相等。

>>> dist, count_edges = np.histogram(count, bins=20, range=(0,20), normed=True)
>>> poisson = stats.poisson.pmf(x, _lambda)
>>> np.max(np.abs(dist-possion)) # 最大偏差很小，符合泊松分布
0.0088356241037075706

还能够换一个角度看随机事件的分布问题。咱们能够观察相邻两个事件之间的时间间隔的分布状况，或者隔k个事件的时间间隔的分布状况。根据几率论，事件之间的时间间隔应符合伽玛分布，因为时间间隔能够是任意数值，所以伽玛分布是一种连续几率分布。伽玛分布的几率密度函数公式以下，它描述第个事件发生所需的等待时间的几率分布。 $\Gamma(k)$ 是伽玛函数，当为整数时，它的值和的阶乘相等。

$f(X;k, \lambda)=\frac{X^{(k-1)}\lambda^{k} e^{(-\lambda X)}}{\Gamma(k)}$

下面的程序模拟了事件的时间间隔的伽玛分布，其结果如【图:模拟伽玛分布】所示。图中的观察时间为1000秒，平均每秒产生10个事件。左图中“k=1”，它表示相邻两个事件间的间隔的分布，而“k=2”则表示相隔一个事件的两个事件间的间隔的分布，能够看出它们都符合伽玛分布。

03-scipy/scipy_gamma_sim.py

模拟伽玛分布

下面咱们直接在IPython中模拟伽玛分布。首先在10000秒以内产生100000个随机事件发生的时刻。所以事件的平均发生次数为每秒10次：

>>> _lambda = 10
>>> time = 10000
>>> t = np.random.rand(_lambda*time)*time

为了计算事件先后的时间间隔，须要先对随机时刻进行排序：

>>> t.sort()

而后分别计算“k=1”和“k=2”时的时间间隔：

>>> s1 = t[1:] - t[:-1]    #相邻两事件的时间间隔
>>> s2 = t[2:] - t[:-2]    #相隔一个事件的两个事件的时间间隔

对s1和s2分别调用histogram()进行几率统计，设置normed为True能够直接统计几率密度：

>>> dist1, x1 = np.histogram(s1, bins=100, normed=True)
>>> dist2, x2 = np.histogram(s2, bins=100, normed=True)

histogram()返回的第二个值为统计区间的边界，下面gamma.pdf()计算伽玛分布的几率密度时，使用各个区间的中值进行计算。pdf()的第二个参数为k值，scale参数为 $1/\lambda$ ：

>>> gamma1 = stats.gamma.pdf((x1[:-1]+x1[1:])/2, 1, scale=1.0/_lambda)
>>> gamma2 = stats.gamma.pdf((x2[:-1]+x2[1:])/2, 2, scale=1.0/_lambda)
>>> np.max(np.abs(gamma1 - dist1))
0.13557317865888141
>>> np.max(np.abs(gamma2 - dist2))
0.087375030861794656

因为几率密度函数的值自己比较大，所以上面的偏差已经很小了：

>>> np.max(gamma1), np.max(gamma2)
(9.3483221580498537, 3.6767953241013656)