统计-stats

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

Footnotes

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

连续和离散概率分布

可以使用下面的语句获得stats模块中所有的连续随机变量:

>>> 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.py

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

下面以正规分布为例,简单地介绍随机变量的用法。下面的语句获得缺省正规分布的随机变量的期望值和方差,我们看到缺省情况下它是一个均值为0,方差为1的随机变量:

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

norm可以像函数一样调用,通过loc和scale参数可以指定随机变量的偏移和缩放参数。对于正态分布的随机变量来说,这两个参数相当于指定其期望值和标准差[2]

>>> 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]))
/tech/static/books/scipy/_images//scipy_stats.png

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

有些随机分布除了loc和scale参数之外,还需要额外的形状参数。例如伽玛分布可用于描述等待k个独立的随机事件发生所需的时间,k就是伽玛分布的形状参数。下面计算形状参数k为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

二项、泊松、伽玛分布

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

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

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

例如可以通过二项分布的概率质量公式计算投掷5次骰子出现3次6点的概率。投掷一次骰子,点数为6的概率(即试验成功的概率)为p=1/6,试验次数为n=5。使用二项分布的概率质量函数pmf()可以很容易计算出现k次6点的概率。和概率密度函数pdf()类似,pmf()的第一个参数为随机变量的取值,后面的参数为描述随机分布所需的参数。对于二项分布来说,参数分别为np,而取值范围则为0到n之间的整数。下面的程序计算k为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%。

在二项分布中,如果试验次数n很大,而每次试验成功的概率p很小,其乘积n p比较适中时,那么试验成功的次数的概率可以用泊松分布近似描述。

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

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

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

03-scipy/scipy_binom_poisson.py

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

程序中事件平均发生率\lambda恒等于10。根据二项分布的试验次数n,计算每次事件出现的概率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
/tech/static/books/scipy/_images//scipy_binom_poisson.png

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

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

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

03-scipy/scipy_poisson_sim.py

模拟泊松分布

/tech/static/books/scipy/_images//scipy_poisson_sim.png

模拟泊松分布

由于上面的程序中包含了许多绘图方面的代码,下面我们直接在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个事件的时间间隔的分布情况。根据概率论,事件之间的时间间隔应符合伽玛分布,由于时间间隔可以是任意数值,因此伽玛分布是一种连续概率分布。伽玛分布的概率密度函数公式如下,它描述第k个事件发生所需的等待时间的概率分布。\Gamma(k)是伽玛函数,当k为整数时,它的值和k的阶乘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

模拟伽玛分布

/tech/static/books/scipy/_images//scipy_gamma_sim.png

模拟伽玛分布

下面我们直接在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)

內容目录

上一个主题

图像处理-ndimage

下一个主题

嵌入C语言程序-weave

本页

loading...