3.4.4 二项分布、泊松分布、伽玛分布
本节用几个实例程序对概率论中的二项分布、泊松分布以及伽玛分布进行一些实验和讨论。
二项分布是最重要的离散概率分布之一。假设有一种只有两个结果的试验，其成功概率为p，那么二项分布描述了进行n次这样的独立试验，成功k次的概率。二项分布的概率质量函数公式如下：
例如，可以通过二项分布的概率质量公式计算投掷5次骰子出现3次6点的概率。投掷一次骰子，点数为6的概率（即试验成功的概率）为p=1/6，试验次数为n=5。使用二项分布的概率质量函数pmf( )可以很容易计算出现k次6点的概率。和概率密度函数pdf( )类似，pmf( )的第一个参数为随机变量的取值，后面的参数为描述随机分布所需的参数。对于二项分布来说，参数分别为n和p，而取值范围则为0到n之间的整数。下面的程序计算k为0到6时对应的概率：
stats.binom.pmf(range(6), 5, 1/6.0)
array([ 4.01877572e-01, 4.01877572e-01, 1.60751029e-01,
3.21502058e-02, 3.21502058e-03, 1.28600823e-04])
由结果可知：出现0或1次6点的概率为40.2%，而出现3次6点的概率为3.215%。
在二项分布中，如果试验次数n很大，而每次试验成功的概率p很小，乘积np比较适中，那么试验成功次数的概率可以用泊松分布近似描述。
在泊松分布中使用λ描述单位时间（或单位面积）中随机事件的平均发生率。如果将二项分布中的试验次数n看作单位时间中所做的试验次数，那么它和事件出现的概率p的乘积就是事件的平均发生率λ，即λ=n·p。泊松分布的概率质量函数公式如下：
下面的程序分别计算二项分布和泊松分布的概率质量函数，结果如图3-14所示。可以看出当n足够大时，二者是十分接近的。程序中的事件平均发生率λ恒等于10。根据二项分布的试验次数n，计算每次事件出现的概率p=λ/n。
图3-14 当n足够大时二项分布和泊松分布近似相等
lambda_ = 10.0
x = np.arange(20)
n1, n2 = 100, 1000
y_binom_n1 = stats.binom.pmf(x, n1, lambda_ / n1)
y_binom_n2 = stats.binom.pmf(x, n2, lambda_ / n2)
y_poisson = stats.poisson.pmf(x, lambda_)
print np.max(np.abs(y_binom_n1 - possion))
print np.max(np.abs(y_binom_n2 - possion))
0.00675531110335
0.000630175404978
泊松分布适合描述单位时间内随机事件发生的次数的分布情况。例如某个设施在一定时间内的使用次数、机器出现故障的次数、自然灾害发生的次数等。
为了加深读者对泊松分布概念的理解，下面我们使用随机数模拟泊松分布，并与概率质量函数进行比较，结果如图3-15所示。图中，每秒内事件的平均发生次数为10，即λ=10。其中下左图的观察时间为1000秒，而下右图的观察时间为50000秒。可以看出观察时间越长，每秒内事件发生的次数越符合泊松分布。
图3-15 模拟泊松分布
np.random.seed(42)
def sim_poisson(lambda_, time):
t = np.random.uniform(0, time, size=lambda_ *
time)❶
count, time_edges = np.histogram(t, bins=time, range=(0, time)) ❷
dist, count_edges = np.histogram(count, bins=20, range=(0, 20), density=True) ❸
x = count_edges[:-1]
poisson = stats.poisson.pmf(x, lambda_)
return x, poisson, dist
lambda_ = 10
times = 1000, 50000
x1, poisson1, dist1 = sim_poisson(lambda_, times[0])
x2, poisson2, dist2 = sim_poisson(lambda_, times[1])
max_error1 = np.max(np.abs(dist1 - poisson1))
max_error2 = np.max(np.abs(dist2 - poisson2))
print "time={}, max_error={}".format(times[0], max_error1)
print "time={}, max_error={}".format(times[1], max_error2)
time=1000, max_error=0.019642302016
time=50000, max_error=0.00179801289496
❶可以用NumPy的随机数生成函数uniform( )产生平均分布于0到time之间的lambda_* time个事件所发生的时刻。❷用histogram( )可以统计数组t中每秒之内的事件发生的次数count，根据泊松分布的定义，count数组中的数值的分布情况应该符合泊松分布。❸接下来统计事件次数在0到20区间内的概率分布。当histogram( )的density参数为True时，结果和概率质量函数相等。
还可以换个角度看随机事件的分布问题。我们可以观察相邻两个事件之间的时间间隔的分布情况，或者隔k个事件的时间间隔的分布情况。根据概率论，事件之间的时间间隔应符合伽玛分布，由于时间间隔可以是任意数值，因此伽玛分布是连续概率分布。伽玛分布的概率密度函数公式如下，它描述第k个事件发生所需的等待时间的概率分布。Γ(k)是伽玛函数，当k为整数时，它的值和k的阶乘k！相等。
下面的程序模拟了事件的时间间隔的伽玛分布，结果如图3-16所示。图中的观察时间为1000秒，平均每秒产生10个事件。下左图中k=1，它表示相邻两个事件间的间隔的分布，而k=2则表示相隔一个事件的两个事件间的间隔的分布，可以看出它们都符合伽玛分布。
图3-16 模拟伽玛分布
def sim_gamma(lambda_, time, k):
t = np.random.uniform(0, time, size=lambda_ *
time) ❶
t.sort( ) ❷
interval = t[k:] - t[:-k] ❸
dist, interval_edges = np.histogram(interval, bins=100, density=True) ❹
x = (interval_edges[1:] + interval_edges[:-1])/2 ❺
gamma = stats.gamma.pdf(x, k, scale=1.0/lambda_) ❺
return x, gamma, dist
lambda_ = 10.0
time = 1000
ks = 1, 2
x1, gamma1, dist1 = sim_gamma(lambda_, time, ks[0])
x2, gamma2, dist2 = sim_gamma(lambda_, time, ks[1])
❶首先在1000秒之内产生10000个随机事件发生的时刻。因此事件的平均发生次数为每秒10次。❷为了计算事件前后的时间间隔，需要先对随机时刻进行排序，❸然后再计算k个事件之间的时间间隔。❹对该时间间隔调用histogram( )进行概率统计，设置density为True可以直接计算概率密度。histogram( )返回的第二个值为统计区间的边界，❺接下来用gamma.pdf( )计算伽玛分布的概率密度时，使用各个区间的中值进行计算。pdf( )的第二个参数为k值，scale参数为1/λ。
接下来我们看一道关于伽玛分布的概率题：有A和B两路公交车，平均发车间隔时间分别是5分钟和10分钟，某乘客在站点S可以任意选择两者之一乘坐，假设A和B到达S的时刻无法确定，计算该乘客的平均等待公交车的时间。
可以将“假设A和B到达S的时刻无法确定”理解为公交车到达S站点的时刻是完全随机的，因此单位时间之内到达S站点的公交车次数符合泊松分布，而前后两辆公交车的时间差符合k=1的伽玛分布。下面我们先用随机数模拟的方法求出近似解，然后推导出解的公式。
T = 100000
A_count = T / 5
B_count = T / 10
A_time = np.random.uniform(0, T, A_count) ❶
B_time = np.random.uniform(0, T, B_count)
bus_time = np.concatenate((A_time, B_time)) ❷
bus_time.sort( )
N = 200000
passenger_time = np.random.uniform(bus_time[0], bus_time[-1], N) ❸
idx = np.searchsorted(bus_time, passenger_time) ❹
np.mean(bus_time[idx] - passenger_time) *
60 ❺
199.12512768644049
模拟的总时间为T分钟，在这段之间之内，应该有A_count次A路公交车和B_count次B路公交车到达S站点。❶可以用均匀分布uniform( )产生两路公交车到达S站点的时刻，❷将这两个保存时刻的数组连接起来，并进行排序。
❸在第一趟和最后一趟公交车的到达时间之间，产生乘客随机到达S站点的时刻数组。❹在已经排序的公交车到达时刻数组bus_time中使用二分法搜索每个乘客到达时刻所在的下标数组idx。❺bus_time[idx]就是乘客到达车站之后第一个到达车站的公交车的时刻，因此只需要计算其差值，并求平均值即可。通过随机数模拟得出的平均等待时间约为200秒。
将A和B两路汽车一起考虑，前后两个车次的平均间隔也为200秒。这似乎有些不可思议，直觉上我们可能期待一个小于平均间隔的等待时间。
np.mean(np.diff(bus_time)) *
60
199.98208112933918
这是因为存在观察者偏差，即会有更多的乘客出现在时间间隔较长的时间段。我们可以想象如果公交车因为事故晚点很长时间，那么通常车站上会挤满等待的人。在图3-17（上）中，蓝色竖线代表公交车的到站时刻，红色竖线代表乘客的到站时刻。可以看出，两条蓝色竖线之间的距离越大，其间的红色竖线就会越多。图3-17（下）的横轴是前后两辆公交车的时间差，纵轴是这段时间差之内的等待人数，可以看出二者成正比关系。
图3-17 观察者偏差
通过以上分析，不难写出计算平均等待时间的计算公式：
在公式中，x是两辆公交车之间的间隔时间，f(x)dx是时间间隔为x出现的概率。由于观察者效应，乘客出现在较长时间间隔的概率也较大，因此xf(x)dx可以看作与乘客出现在时间间隔为x时段的概率成比例的量，分母的积分将其归一化。而分子中的x/2是在该时间间隔段到达车站所需的平均等待时间。下面我们计算该公式，由图3-17可知，公交车的间隔几乎不会超过30分钟，因此虽然公式中的积分上限为+∞，但在实际计算时只需要指定一个较大的数即可。在本章后续的小节中会详细介绍数值积分quad( )的用法。
from scipy import integrate
t = 10.0 / 3 # 两辆公交车之间的平均时间间隔
bus_interval = stats.gamma(1, scale=t)
n, _ = integrate.quad(lambda x: 0.5 *
x *
x *
bus_interval.pdf(x), 0, 1000)
d, _ = integrate.quad(lambda x: x *
bus_interval.pdf(x), 0, 1000)
n / d *
60
200.0