SciPy的stats模块包含了多种概率分布的随机变量[1],随机变量分为连续和离散两种。所有的连续随机变量都是rv_continuous的派生类的对象,而所有的离散随机变量都是rv_discrete的派生类的对象。
Footnotes
[1] | 本节中的随机变量是指概率论中的概念,不是Python中的变量 |
连续和离散概率分布
可以使用下面的语句获得stats模块中所有的连续随机变量:
连续随机变量对象都有如下方法:
- rvs:对随机变量进行随机取值,可以通过size参数指定输出的数组的大小。
- pdf:随机变量的概率密度函数。
- cdf:随机变量的累积分布函数,它是概率密度函数的积分。
- sf:随机变量的生存函数,它的值是1-cdf(t)。
- ppf:累积分布函数的反函数。
- stat:计算随机变量的期望值和方差。
- fit:对一组随机取样进行拟合,找出最适合取样数据的概率密度函数的系数。
概率密度函数、直方图统计和累积分布函数
下面以正规分布为例,简单地介绍随机变量的用法。下面的语句获得缺省正规分布的随机变量的期望值和方差,我们看到缺省情况下它是一个均值为0,方差为1的随机变量:
norm可以像函数一样调用,通过loc和scale参数可以指定随机变量的偏移和缩放参数。对于正态分布的随机变量来说,这两个参数相当于指定其期望值和标准差[2]:
下面调用随机变量X的rvs()方法,得到包含一万次随机取样值的数组x,然后调用NumPy的mean()和var()计算此数组的均值和方差,其结果符合随机变量X的特性:
我们也可以使用fit()方法对随机取样序列x进行拟合,它返回的是与随机取样值最吻合的随机变量的参数:
接下来比较随机变量X的概率密度函数和对数组x进行直方图统计的结果:
其中histogram()对数组x进行直方图统计,它将数组x的取值范围分为100个区间,并统计x中的每个值落入各个区间中的次数。histogram()返回两个数组p和t2,其中p表示各个区间的取样值出现的频数,由于normed参数为True,因此p的值是正规化之后的结果。t2表示区间,由于其中包括区间起点和终点,因此t2的长度为101。【图:正规分布的概率密度函数(左)和累积分布函数(右)】(左)显示了概率密度函数和直方图统计的结果,可以看出二者是一致的。
下面的程序绘制随机变量X的累积分布函数和数组p的累加结果,其结果如【图:正规分布的概率密度函数(左)和累积分布函数(右)】(右)所示。
有些随机分布除了loc和scale参数之外,还需要额外的形状参数。例如伽玛分布可用于描述等待个独立的随机事件发生所需的时间,就是伽玛分布的形状参数。下面计算形状参数为1和2时的伽玛分布的期望值和方差:
伽玛分布的尺度参数和随机事件发生的频率相关,由scale参数指定:
根据伽玛分布的数学定义可知其期望值为,方差为。上面的程序验证了这两个公式。
当随机分布有额外的形状参数时,它所对应的rvs()、pdf()等方法都会增加额外的参数接收形状参数。例如下面的程序调用rvs()对的伽玛分布取4个随机值:
接下来调用pdf()查看上面4个随机值所对应的概率密度:
也可以先创建将形状参数和尺度参数固定的随机变量,然后再调用其pdf()计算概率密度:
当分布函数的值域为离散时我们称之为离散概率分布。例如投掷有六个面的骰子时,只能获得1到6的整数,因此所得到的概率分布为离散的。对于离散随机分布,通常使用概率质量函数(PMF)描述其分布情况。
在stats库中所有描述离散分布的随机变量都从rv_discrete类继承。也可以直接用rv_discrete类自定义离散概率分布。例如假设有一个不均匀的骰子,它的各点出现的概率不相等。我们可以用下面的数组x保存骰子的所有可能值,数组p保存每个值出现的概率:
于是可以用下面的语句定义表示这个特殊骰子的随机变量,并调用其rvs()方法投掷此骰子20次,获得符合概率p的随机数:
Footnotes
[2] | 标准差是方差的算术平方根,因此标准差为2.0时,方差为4.0 |
二项、泊松、伽玛分布
本节用几个实例程序对概率论中的二项分布、泊松分布以及伽玛分布进行一些实验和讨论。
二项分布是最重要的离散概率分布之一。假设有一种只有两个结果的试验,其成功概率为,那么二项分布描述了进行次这样的独立试验,成功次的概率。二项分布的概率质量函数公式如下:
例如可以通过二项分布的概率质量公式计算投掷5次骰子出现3次6点的概率。投掷一次骰子,点数为6的概率(即试验成功的概率)为,试验次数为。使用二项分布的概率质量函数pmf()可以很容易计算出现次6点的概率。和概率密度函数pdf()类似,pmf()的第一个参数为随机变量的取值,后面的参数为描述随机分布所需的参数。对于二项分布来说,参数分别为和,而取值范围则为0到之间的整数。下面的程序计算为0到6所对应的概率:
由结果可知:出现0或1次6点的概率为40.2%,而出现3次6点的概率为3.215%。
在二项分布中,如果试验次数很大,而每次试验成功的概率很小,其乘积比较适中时,那么试验成功的次数的概率可以用泊松分布近似描述。
在泊松分布中使用描述单位时间(或单位面积)中随机事件的平均发生率。如果我们将二项分布中的试验次数看作单位时间中所做的试验次数,那么它和事件出现的概率的乘积就是事件的平均发生率,即。泊松分布的概率质量函数公式如下:
下面的程序分别计算二项分布和泊松分布的概率质量函数,其结果如【图:当n足够大时二项分布和泊松分布近似相等】所示,可以看出当足够大时,二者是十分接近的。
03-scipy/scipy_binom_poisson.py
比较二项分布和泊松分布的概率质量函数
程序中事件平均发生率恒等于10。根据二项分布的试验次数,计算每次事件出现的概率。程序中运算部分大致如下:
泊松分布适合描述单位时间内随机事件发生的次数的分布情况。例如某一设施在一定时间内的使用次数,机器出现故障的次数,自然灾害发生的次数等等。
为了加深读者对泊松分布概念的理解,下面我们使用随机数模拟泊松分布,并与其概率质量函数进行比较,其结果如【图:模拟泊松分布】所示。图中,每秒内事件的平均发生次数为10,即。其中左图的观察时间为1000秒,而右图的观察时间为50000秒。可以看出观察时间越长,每秒内事件发生的次数越符合泊松分布。
由于上面的程序中包含了许多绘图方面的代码,下面我们直接在IPython中介绍泊松分布的模拟过程。首先定义事件发生率和观察时间:
可以用NumPy的随机数生成函数rand(),产生平均分布于0到time之间的_lambda*time个事件所发生的时刻。由于rand()产生的是0到1之间的平均分布的随机数,因此需要对其结果扩大time倍:
用histogram()可以统计数组t中每秒之内的事件发生的次数count:
根据泊松分布的定义,count数组中的数值的分布情况应该符合泊松分布。下面统计事件次数在0到20区间内的概率分布。当histogram()的normed参数为True并且每个统计区间的长度为1时,其结果和概率质量函数相等。
还可以换一个角度看随机事件的分布问题。我们可以观察相邻两个事件之间的时间间隔的分布情况,或者隔k个事件的时间间隔的分布情况。根据概率论,事件之间的时间间隔应符合伽玛分布,由于时间间隔可以是任意数值,因此伽玛分布是一种连续概率分布。伽玛分布的概率密度函数公式如下,它描述第个事件发生所需的等待时间的概率分布。是伽玛函数,当为整数时,它的值和的阶乘相等。
下面的程序模拟了事件的时间间隔的伽玛分布,其结果如【图:模拟伽玛分布】所示。图中的观察时间为1000秒,平均每秒产生10个事件。左图中“k=1”,它表示相邻两个事件间的间隔的分布,而“k=2”则表示相隔一个事件的两个事件间的间隔的分布,可以看出它们都符合伽玛分布。
模拟伽玛分布
下面我们直接在IPython中模拟伽玛分布。首先在10000秒之内产生100000个随机事件发生的时刻。因此事件的平均发生次数为每秒10次:
为了计算事件前后的时间间隔,需要先对随机时刻进行排序:
然后分别计算“k=1”和“k=2”时的时间间隔:
对s1和s2分别调用histogram()进行概率统计,设置normed为True可以直接统计概率密度:
histogram()返回的第二个值为统计区间的边界,下面gamma.pdf()计算伽玛分布的概率密度时,使用各个区间的中值进行计算。pdf()的第二个参数为k值,scale参数为:
由于概率密度函数的值本身比较大,因此上面的误差已经很小了: